#This program solves the Navier Stokes equation. It uses the method of
#constructed solutions to test for convergence, which should be 2nd order 

#from __future__ import division
import numpy as np
import numpy.fft as fft
from numpy import pi 
from pyfft.cl import Plan
import pyopencl as cl
import pyopencl.array as cl_array
from pyopencl.elementwise import ElementwiseKernel
import time
import sys
#import matplotlib.delaunay as triang
import meshpy.triangle as triangle
from pyvisfile.vtk import write_structured_grid
from pyvisfile.vtk import (
           UnstructuredGrid, DataArray,
           AppendedDataXMLGenerator, VTK_TRIANGLE,
           VTK_LINE, VTK_VERTEX, VF_LIST_OF_VECTORS, VF_LIST_OF_COMPONENTS)
#np.set_printoptions(threshold=sys.maxint)
from warnings import filterwarnings
filterwarnings("ignore", module="pyopencl.cache", lineno=336)
filterwarnings("ignore", module="pyopencl", lineno=163)
TIMING = 0 
WGPSIZE_1 = 16
WGPSIZE_2 = 16 
b_WGPSIZE = 4    
PR_INT = 400
DISP_INT = 400
#PR_INT = 800
#DISP_INT = 800
PR_INT = 4000
DISP_INT = 4000
#PR_INT = 8000
#DISP_INT = 8000
SAVE = 1
start_size = 9
end_size = 9

#f = open('/home/tfai/comfile', 'w')

affine_part = np.array([3,3,1,1], dtype=np.float64)/8.
r_vec = np.array([-2,-2,2,2], dtype=np.float64)/8.
q_vec = np.array([-1,1,1,-1], dtype=np.float64)/8.

def myq(r):
    return np.sqrt(1+4*r*(1-r))

def spreadF(F_1, F_2, X_1, X_2, b_q, dx, dtheta): #spread Lagrangian forces
    
    N = np.int(np.round(1/dx,0))
    c = dtheta/(dx*dx);
    f_1 = np.zeros((N,N))
    f_2 = np.zeros((N,N))

    for k in range(0, b_q):
        s_1_stag = (X_1[k]+0.5)/dx
        s_2_stag = (X_2[k]+0.5)/dx        
        i_1_stag = np.int(np.floor(s_1_stag))
        i_2_stag = np.int(np.floor(s_2_stag))
        r_1_stag = s_1_stag-i_1_stag
        r_2_stag = s_2_stag-i_2_stag
        s_1_c = (X_1[k]+0.5-dx/2)/dx
        s_2_c = (X_2[k]+0.5-dx/2)/dx        
        i_1_c = np.int(np.floor(s_1_c))
        i_2_c = np.int(np.floor(s_2_c))
        r_1_c = s_1_c-i_1_c
        r_2_c = s_2_c-i_2_c

	w_1 = (
		(affine_part + r_1_stag*r_vec + myq(r_1_stag)*q_vec).reshape(1,-1)
		* (affine_part + r_2_c*r_vec + myq(r_2_c)*q_vec)[:,np.newaxis])
	w_2 = (
		(affine_part + r_1_c*r_vec + myq(r_1_c)*q_vec).reshape(1,-1)
		* (affine_part + r_2_stag*r_vec + myq(r_2_stag)*q_vec)[:,np.newaxis])

	if False:
	    ind_1_lr = tile((fmod(range(i_1_stag-1,i_1_stag+3),N)),(4,1))
	    ind_2_lr = tile(fmod(range(i_2_c-1,i_2_c+3),N), (4,1)).T
	    ind_1_tb = tile((fmod(range(i_1_c-1,i_1_c+3),N)),(4,1))
	    ind_2_tb = tile(fmod(range(i_2_stag-1,i_2_stag+3),N), (4,1)).T

	    f_1[ind_2_lr,ind_1_lr] += c*F_1[k]*w_1
	    f_2[ind_2_tb,ind_1_tb] += c*F_2[k]*w_2
	else:

	    ind_1_lr = (np.arange(i_1_stag-1,i_1_stag+3, dtype=np.int32) % N)[np.newaxis,:]
	    ind_2_lr = (np.arange(i_2_c-1,i_2_c+3, dtype=np.int32) % N)[:,np.newaxis]
	    ind_1_tb = (np.arange(i_1_c-1,i_1_c+3, dtype=np.int32) % N)[np.newaxis,:]
	    ind_2_tb = (np.arange(i_2_stag-1,i_2_stag+3, dtype=np.int32) % N)[:,np.newaxis]

	    f_1[ind_2_lr,ind_1_lr] += c*F_1[k]*w_1
	    f_2[ind_2_tb,ind_1_tb] += c*F_2[k]*w_2

    return f_1, f_2

def retriang(flag):
    global b_q, b_add, N_edge, N_tri, TV
    global TV_1_gpu, TV_2_gpu, TV_3_gpu
    global X_1_gpu, X_2_gpu, X_1_tilde_gpu, X_2_tilde_gpu
    global F_1_gpu, F_2_gpu, U_1_gpu, U_2_gpu
    global Rho_gpu, Mu_gpu
    global b_WGPSIZE_2, b_WGPSIZE_3
    global X_arr_1_gpu, X_arr_2_gpu
    global tmpRho_gpu, tmpMu_gpu
    global edges, edges_array
    #using meshpy
    X_arr = np.hstack((X_1_gpu.get()[0:b_q],X_2_gpu.get()[0:b_q])).copy()
    X_arr = np.vstack((X_arr,np.hstack((X_1_gpu.get()[b_q+N_teth:b_q+N_teth+b_add],X_2_gpu.get()[b_q+N_teth:b_q+N_teth+b_add]))))
    #X_arr = np.vstack((X_arr,np.hstack((X_1_gpu.get()[b_q+N_teth:b_q+N_teth+N_edge],X_2_gpu.get()[b_q+N_teth:b_q+N_teth+N_edge]))))
    info = triangle.MeshInfo()
    #info.set_points(list(X_arr))
    #info.set_points(list(X_arr), range(2,b_q+2))
    info.set_points(list(X_arr), range(2,b_q+2)+range(2+b_q+N_teth,2+b_q+N_teth+b_add))
    #info.set_points(list(X_arr), range(2,b_q+N_edge+2))
    #L1 = round_trip_connect(0, N_s-1)
    #L2 = round_trip_connect(N_s*(N_r-1),N_s*N_r-1)
    #info.set_facets(L1+L2)
    info.set_facets(edges)
    #mesh = triangle.build(info, allow_volume_steiner=False,
    #   quality_meshing=False, allow_boundary_steiner=False)
    #mesh = triangle.build(info,
    #   min_angle=2.0, allow_boundary_steiner=False)

    #print TV.shape
    #print np.asarray(edges).shape
    #print X_arr.shape

    #mesh = triangle.build(info,
    #   min_angle=1.0, max_volume=20.0*dx*dx)
    #mesh = triangle.build(info,
    #   min_angle=1.0)
    mesh = triangle.build(info,
       min_angle=1.0, max_volume=10.0*dx*dx)
    temp = len(mesh.points)
    markers = np.array(mesh.point_markers)
    points = np.array(mesh.points)
    #mesh.points = np.vstack((mesh.points[markers != 1],mesh.points[markers == 1]))

    TV = np.array(mesh.elements).copy()

    ##try to get edges manually
    #edge_to_el = {}
    #for el_nr, el in enumerate(TV):
    #    i, j, k = el

    #    for e in [
    #            frozenset([i, j]),
    #            frozenset([j, k]),
    #            frozenset([i, k]),
    #            #[i, j],
    #            #[j, k],
    #            #[i, k],
    #            ]:
    #        edge_to_el.setdefault(e, []).append(el_nr)

    #edges = []
    #for edge, next_door_elements in edge_to_el.iteritems():
    #    if len(next_door_elements) == 1:
    #        #print edge
    #        #edges.append(edge)
    #        edges.append(list(edge))

    #print edges
    #edges = np.asarray(edges)
    myedges = np.asarray(mesh.facets)
    #edges = np.sort(edges)
    #edges2 = np.sort(edges2)
    #sind = np.lexsort((edges[:,1],edges[:,0]))
    #sind2 = np.lexsort((edges2[:,1],edges2[:,0]))
    ##print np.max(np.abs(np.sort(edges)-np.sort(edges2))) 
    #print np.max(np.abs(edges[sind,:]-edges2[sind2,:]))
    ##print np.sort(edges)
    ##print np.sort(edges2)
    #print edges[sind,:]
    #print edges2[sind2,:]
    #sys.exit()
    inds = np.nonzero(markers == 1)[0] 
    print 'number of ones is', len(inds)
    #print 'locations of ones are', inds
    inds = np.nonzero(markers == 0)[0] 
    print 'number of zeros is', len(inds)
    #print 'locations of zeros are', inds
    print 'number of vertices is', temp
    print 'number of original vertices is', b_q
    #print edges
    #for i in range(0,len(inds)):
    #    print inds[i]
    #    print points[inds[i],:]
    #    print edges[edges[:,0] == inds[i],:]
    #    print edges[edges[:,1] == inds[i],:]
    #    points = np.vstack((points[:inds[i],:],points[inds[i]+1:,:],points[inds[i],:]))
    #    #print (markers[:inds[i]][:,np.newaxis]).shape
    #    #print (markers[inds[i]:][:,np.newaxis]).shape
    #    #print points[inds[i]].shape
    #    markers = np.vstack((markers[:inds[i]][:,np.newaxis],markers[inds[i]+1:][:,np.newaxis],markers[inds[i]])).flatten(1)
    #    TV_old = TV.copy()
    #    TV[TV_old > inds[i]] = TV[TV_old > inds[i]]-1 
    #    TV[TV_old == inds[i]] = temp-1 
    #    edges_old = edges.copy() 
    #    edges[edges_old > inds[i]] = edges[edges_old > inds[i]]-1 
    #    edges[edges_old == inds[i]] = temp-1 
    #    for j in range(i+1,len(inds)):
    #        if (inds[j] > inds[i]): 
    #            inds[j] = inds[j] - 1
    ##print edges

    #number_ones = 0
    #ind_ones = 0
    #while (ind_ones >= 0):
    #    if (markers[ind_ones] == 1):
    #        number_ones += 1
    #    else:
    #        ind_ones = -1

    #print markers.shape
    #print (markers[markers != 1]).shape
    #print (markers[markers == 1]).shape

    inds = np.nonzero(markers <= 1)[0] 

    #print inds
    print len((np.nonzero(markers > 1)[0]))
    #print markers
    #print b_q
    #print N_teth
    #print np.asarray(mesh.facets)
    #print np.asarray(mesh.elements)
    #sys.exit()

    #b_add = np.array(mesh.points).shape[0]-b_q
    #if (b_add != 0):
    if (len(inds) != 0):
        print 'flag is', flag
        print 'new point added'
        b_add_old = b_add
        #b_add = len(inds)
        b_add = b_add_old+len(inds)
        X_arr = np.hstack((X_1_gpu.get(),X_2_gpu.get())).copy()
        X_arr = np.vstack((X_arr, points[inds,:])).copy()
        if (flag == 0):
            X_arr = np.hstack((X_1_tilde_gpu.get(),X_2_tilde_gpu.get())).copy()
            X_arr = np.vstack((X_arr, points[inds,:])).copy()
        X_arr_1_gpu = cl_array.to_device(queue, X_arr[:,0][:,np.newaxis].astype(np.float64))
        X_arr_2_gpu = cl_array.to_device(queue, X_arr[:,1][:,np.newaxis].astype(np.float64))
        #tmpMu = np.zeros((b_add,1))
        #tmpRho = np.zeros((b_add,1))
        tmpMu = np.zeros((b_add-b_add_old,1))
        tmpRho = np.zeros((b_add-b_add_old,1))
        myRho = Rho_gpu.get()
        myMu = Mu_gpu.get()
        tmpMu_gpu = cl_array.to_device(queue, tmpMu.astype(np.float64))
        tmpRho_gpu = cl_array.to_device(queue, tmpRho.astype(np.float64))
        #calc_coefs_retriang_kernel(queue, (N_tri, 1), (b_WGPSIZE_3,1), X_arr_1_gpu.data, X_arr_2_gpu.data, TV_1_gpu.data, TV_2_gpu.data, TV_3_gpu.data, Rho_gpu.data, Mu_gpu.data, b_add, b_q+N_teth+b_add_old, tmpRho_gpu.data, tmpMu_gpu.data) 
        calc_coefs_retriang_kernel(queue, (N_tri, 1), (b_WGPSIZE_3,1), X_arr_1_gpu.data, X_arr_2_gpu.data, TV_1_gpu.data, TV_2_gpu.data, TV_3_gpu.data, Rho_gpu.data, Mu_gpu.data, b_add-b_add_old, b_q+N_teth+b_add_old, tmpRho_gpu.data, tmpMu_gpu.data) 
        tmpMu = tmpMu_gpu.get()
        tmpRho = tmpRho_gpu.get()
        #myMu = np.vstack((Mu, tmpMu)).copy() 
        #myRho = np.vstack((Rho, tmpRho)).copy() 
        myMu = np.vstack((myMu, tmpMu)).copy() 
        myRho = np.vstack((myRho, tmpRho)).copy() 
        print 'zeros is tmpRho is', len(np.nonzero(tmpRho == 0.0)[0])
        Rho_gpu = cl_array.to_device(queue, myRho.astype(np.float64)) 
        Mu_gpu = cl_array.to_device(queue, myMu.astype(np.float64)) 
        #X_arr = np.hstack((X_1_gpu.get(),X_2_gpu.get())).copy()
        #X_arr = np.vstack((X_arr[0:b_q+N_teth,:], points[inds,:])).copy()
        X_arr = np.hstack((X_1_gpu.get(),X_2_gpu.get())).copy()
        X_arr = np.vstack((X_arr, points[inds,:])).copy()
        X_1_gpu = cl_array.to_device(queue, X_arr[:,0][:,np.newaxis].astype(np.float64))
        X_2_gpu = cl_array.to_device(queue, X_arr[:,1][:,np.newaxis].astype(np.float64))
        X_1_tilde_gpu = cl_array.empty_like(X_1_gpu)
        X_2_tilde_gpu = cl_array.empty_like(X_1_gpu)
        U_1_gpu = cl_array.empty_like(X_1_gpu) 
        U_2_gpu = cl_array.empty_like(X_1_gpu) 
        if (np.mod(b_q+N_teth+b_add, 8) == 0):
            b_WGPSIZE_2 = 8
        elif (np.mod(b_q+N_teth+b_add, 4) == 0):
            b_WGPSIZE_2 = 4
        elif (np.mod(b_q+N_teth+b_add, 2) == 0):
            b_WGPSIZE_2 = 2
        else:
            b_WGPSIZE_2 = 1
        interp_kernel(queue, (b_q+N_teth+b_add, 1), (b_WGPSIZE_2, 1), u_1_gpu.data, u_2_gpu.data, X_1_gpu.data, X_2_gpu.data, N, dx, U_1_gpu.data, U_2_gpu.data)
        add_gpu(X_1_gpu, U_1_gpu, dt, X_1_tilde_gpu)
        add_gpu(X_2_gpu, U_2_gpu, dt, X_2_tilde_gpu)
        F_1_gpu = cl_array.empty_like(X_1_gpu) 
        F_2_gpu = cl_array.empty_like(X_1_gpu) 
        for i in range(0,len(inds)):
            #markers[inds[i]] = N_teth+b_q+i+2
            markers[inds[i]] = N_teth+b_q+b_add_old+i+2
    inds = np.nonzero(markers != np.arange(2,temp+2))[0]
    #TV = np.array(mesh.elements).copy()
    #TV[TV >= b_q] = TV[TV >= b_q] + N_teth
    TV_old = TV.copy()
    #print edges

    #print np.asarray(mesh.facets).shape

    #print edges
    #edges = np.asarray(mesh.facets)
    #edges[edges >= b_q] = edges[edges >= b_q] + N_teth
    #print edges

    #print N_teth
    #print b_q
    #print b_add

    #print edges
    edges_old = myedges.copy()
    N_tri = len(mesh.elements) 
    if (np.mod(N_tri, 16) == 0):
        b_WGPSIZE_3 = 16
    elif (np.mod(N_tri, 8) == 0):
        b_WGPSIZE_3 = 8
    elif (np.mod(N_tri, 4) == 0):
        b_WGPSIZE_3 = 4
    elif (np.mod(N_tri, 2) == 0):
        b_WGPSIZE_3 = 2
    else:
        b_WGPSIZE_3 = 1
    print 'new number of triangles is', N_tri
    for i in range(len(inds)):
        TV[TV_old == inds[i]] = markers[inds[i]]-2 
        myedges[edges_old == inds[i]] = markers[inds[i]]-2 
    #print TV.shape
    #print edges.shape
    #myedges = []
    #edges_array = np.asarray(myedges)
    edges_array = myedges.copy()
    myedges[(myedges > b_q)] = myedges[(myedges > b_q)] - N_teth
    edges = []
    for i in range(0,myedges.shape[0]):
        edges.append((int(myedges[i][0]),int(myedges[i][1]))) 
    #for i in range(0,edges.shape[0]):
    #  if (edges[i,0] >= b_q):
    #      e1 = edges[i,0] - N_teth
    #  else:
    #      e1 = edges[i,0]
    #  if (edges[i,1] >= b_q):
    #      e2 = edges[i,1] - N_teth
    #  else:
    #      e2 = edges[i,1]
    #  myedges.append((int(e1), int(e2)))
    #edges = myedges

    #print np.asarray(edges).shape

    #print edges
    #sys.exit()
    TV_1_gpu = cl_array.to_device(queue, TV[:,0].astype(np.uint32))
    TV_2_gpu = cl_array.to_device(queue, TV[:,1].astype(np.uint32))
    TV_3_gpu = cl_array.to_device(queue, TV[:,2].astype(np.uint32))

def round_trip_connect(start, end):
  result = []
  for i in range(start, end):
    result.append((i, i+1))
  result.append((end, start))
  return result

def spread_routine(q, firstn1, nextn1, firstn2, nextn2, F_1, F_2, X_1, X_2, N, dx, dtheta, f_1, f_2):
    for shift0 in range(0,4):
        for shift1 in range(0,4):
            spread_kernel(q, (N/4, N/4), (WGPSIZE_1, WGPSIZE_2), firstn1.data, nextn1.data, firstn2.data, nextn2.data, F_1.data, F_2.data, X_1.data, X_2.data, N, dx, dtheta, shift0, shift1, f_1.data, f_2.data)

def update_routine(q,  firstn1, nextn1, firstn2, nextn2, X_1, X_2, N, dx):
    for shift0 in range(0,4):
        for shift1 in range(0,4):
            updaten_kernel(q, (N/4, N/4), (WGPSIZE_1, WGPSIZE_2), firstn1.data, nextn1.data, firstn2.data, nextn2.data, X_1.data, X_2.data, N, dx, shift0, shift1)

def calc_transf(dx):
    N = np.int(np.round(1/dx,0))
    k = np.arange(0, N/2+1)
    k = np.hstack([k, np.arange(-N/2+1, 0)])
    k1 = np.tile(k, (N, 1))
    k2 = k1.T.copy()
    L_tf = -4./dx**2 * (np.sin(pi*k1/N)**2+np.sin(pi*k2/N)**2) 
    L_tf[0,0] = 1. #avoid dividing by zero
    return L_tf

def calc_inds(N): #calculate shifted indices
    a1 = range(0,N)
    i_1 = np.tile(a1, (N, 1))
    i_2 = i_1.T.copy()
    a2 = np.hstack([range(1, N), 0])
    ip_1 = np.tile(a2, (N,1))
    ip_2 = ip_1.T.copy()
    a3 = np.hstack([N-1, range(0, N-1)])
    im_1 = np.tile(a3, (N,1))
    im_2 = im_1.T.copy()
    return i_1, i_2, ip_1, ip_2, im_1, im_2

def solve(t_step):
    global dx, dt, mu_0, rho_0, b_q, dtheta, ds, ks, kb, rho_av, N, N_tri
    N = np.int(np.round(1/dx,0))
    copy_gpu(p_gpu, p_last_gpu)

    zero_gpu(en_gpu)
    check_tri_kernel(queue, (N_tri, 1), (b_WGPSIZE_3, 1), X_1_gpu.data, X_2_gpu.data, TV_1_gpu.data, TV_2_gpu.data, TV_3_gpu.data, en_gpu.data) 
    if (en_gpu.get()[0] == 1.0):
        #print 'tangled in beginning'
        #sys.exit()
        retriang(0)

    #div_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), u_1_gpu.data, u_2_gpu.data, dx, N, p_gpu.data)
    #print 'div max is', p_gpu.get().max()

    #print 'u_1 is', u_1_gpu.get().max()
    #print 'u_2 is', u_2_gpu.get().max()
    imag_gpu(u_1_gpu, u_1_tf_gpu)
    imag_gpu(u_2_gpu, u_2_tf_gpu)

    plan.execute(u_1_tf_gpu.data)
    plan.execute(u_2_tf_gpu.data)
    
    interp_kernel(queue, (b_q+N_teth+b_add, 1), (b_WGPSIZE_2, 1), u_1_gpu.data, u_2_gpu.data, X_1_gpu.data, X_2_gpu.data, N, dx, U_1_gpu.data, U_2_gpu.data)
    #interp_kernel(queue, (b_q, 1), (b_WGPSIZE_2, 1), u_1_gpu.data, u_2_gpu.data, X_1_gpu.data, X_2_gpu.data, N, dx, U_1_gpu.data, U_2_gpu.data)
    add_gpu(X_1_gpu, U_1_gpu, dt, X_1_tilde_gpu)
    add_gpu(X_2_gpu, U_2_gpu, dt, X_2_tilde_gpu)

    zero_gpu(F_1_gpu)
    zero_gpu(F_2_gpu)

    zero_gpu(f_1_gpu)
    zero_gpu(f_2_gpu)
    zero_gpu(u_1_tilde_gpu)
    zero_gpu(u_2_tilde_gpu)
    zero_gpu(mu_1_gpu)
    zero_gpu(mu_2_gpu)
    zero_gpu(rho_1_gpu)
    zero_gpu(rho_2_gpu)

    #calc_f_gpu(dt*(t_step-1.0), dx, xx_1_stag_gpu, xx_2_stag_gpu, xx_1_c_gpu, xx_2_c_gpu, f_1_gpu, f_2_gpu) 

    calc_lagF_teth_kernel(queue, (N_teth, 1), (b_WGPSIZE_4, 1), X_1_tilde_gpu.data, X_2_tilde_gpu.data, X_1_teth_gpu.data, X_2_teth_gpu.data, N_pre_teth, k_teth, dx_teth, F_1_gpu.data, F_2_gpu.data)

    update_routine(queue, firstn1_gpu, nextn1_gpu, firstn2_gpu, nextn2_gpu, X_1_tilde_gpu, X_2_tilde_gpu, N, dx)
    spread_routine(queue, firstn1_gpu, nextn1_gpu, firstn2_gpu, nextn2_gpu, F_1_gpu, F_2_gpu, X_1_tilde_gpu, X_2_tilde_gpu, N, dx, dtheta, u_1_tilde_gpu, u_2_tilde_gpu)
    #spread_routine(queue, firstn1_gpu, nextn1_gpu, firstn2_gpu, nextn2_gpu, F_1_gpu, F_2_gpu, X_1_gpu, X_2_gpu, N, dx, dtheta, u_1_tilde_gpu, u_2_tilde_gpu)

    add_gpu(u_2_tilde_gpu, region_gpu, body_f*np.cos(wave_omega*dt*t_step), u_2_tilde_gpu)
    #add_gpu(u_2_tilde_gpu, region_gpu, body_f*np.cos(0.71*dt*t_step), u_2_tilde_gpu)
    #calculate mu, rho on grid

    zero_gpu(en_gpu)
    check_tri_kernel(queue, (N_tri, 1), (b_WGPSIZE_3, 1), X_1_tilde_gpu.data, X_2_tilde_gpu.data, TV_1_gpu.data, TV_2_gpu.data, TV_3_gpu.data, en_gpu.data) 
    if (en_gpu.get()[0] == 1.0):
        retriang(1)
        zero_gpu(en_gpu)
        check_tri_kernel(queue, (N_tri, 1), (b_WGPSIZE_3, 1), X_1_gpu.data, X_2_gpu.data, TV_1_gpu.data, TV_2_gpu.data, TV_3_gpu.data, en_gpu.data) 
        if (en_gpu.get()[0] == 1.0):
            print 'wrong assumption'
            sys.exit()
         
    calc_coefs_kernel(queue, (N_tri, 1), (b_WGPSIZE_3,1), X_1_gpu.data, X_2_gpu.data, TV_1_gpu.data, TV_2_gpu.data, TV_3_gpu.data, Rho_gpu.data, Mu_gpu.data, dx, N, rho_1_gpu.data, rho_2_gpu.data, mu_1_gpu.data, mu_2_gpu.data) 
    #calc_coefs_kernel(queue, (1, 1), (1,1), X_1_gpu.data, X_2_gpu.data, TV_1_gpu.data, TV_2_gpu.data, TV_3_gpu.data, Rho_gpu.data, Mu_gpu.data, dx, N, rho_1_gpu.data, rho_2_gpu.data, mu_1_gpu.data, mu_2_gpu.data) 

    #if SAVE:
    #    copy_gpu(rho_1_gpu, rho_1_save_gpu)
    #    copy_gpu(rho_2_gpu, rho_2_save_gpu)
    #    add_scalar_gpu(rho_1_save_gpu, rho_0, rho_1_save_gpu)
    #    add_scalar_gpu(rho_2_save_gpu, rho_0, rho_2_save_gpu)

    zero_gpu(en_sum_gpu)
    calc_rho_kernel(queue, (N_tri, 1), (b_WGPSIZE_3, 1), X_1_gpu.data, X_2_gpu.data, TV_1_gpu.data, TV_2_gpu.data, TV_3_gpu.data, Rho_gpu.data, en_sum_gpu.data)

    rho_av = (en_sum_gpu.get().sum()+dx*dx*N*N)/(dx*dx*N*N)
    rho_av = (rho_0*N*N+rho_2_gpu.get().sum())/(N*N) 
    #print 'rho_av is', rho_av
    #rho_av = 0
    #rho_av = np.mean(myrho_0*np.exp(-2*pi*xx_2_stag/C_h))

    #add_gpu(u_2_tilde_gpu, rho_0-rho_av+rho_2_gpu, -grav, u_2_tilde_gpu)
    copy_gpu(rho_2_gpu, rho_2_extra_gpu)
    add_scalar_gpu(rho_2_extra_gpu, rho_0-rho_av, rho_2_extra_gpu)
    add_gpu(u_2_tilde_gpu, rho_2_extra_gpu, -grav, u_2_tilde_gpu)

    #print rho_2_gpu.get().max()
    #print mu_1_gpu.get().max()
    #print mu_2_gpu.get().max()
    
    #print u_1_tilde_gpu.get()
    #print u_2_tilde_gpu.get()
    #print 'av force is', np.mean(u_2_tilde_gpu.get())

    #print 'density'
    #print rho_0-rho_av+np.mean(rho_2_gpu.get()) 
    #print (np.sum(Rho)*dtheta+dx**2*rho_0*N**2)/(dx**2*N**2)
    #print mu_1_gpu.get().max()
    #print rho_1_gpu.get().max()
    add_scalar_gpu(rho_1_gpu, rho_0-rho_hat, rho_1_gpu)
    add_scalar_gpu(rho_2_gpu, rho_0-rho_hat, rho_2_gpu)
    #add_scalar_gpu(mu_1_gpu, -mu_bar, mu_1_gpu)
    #add_scalar_gpu(mu_2_gpu, -mu_bar, mu_2_gpu)
    #add_scalar_gpu(mu_1_gpu, mu_0, mu_1_gpu)
    #add_scalar_gpu(mu_2_gpu, mu_0, mu_2_gpu)
    #add_scalar_gpu(rho_1_gpu, -rho_bar, rho_1_gpu)
    #add_scalar_gpu(rho_2_gpu, -rho_bar, rho_2_gpu)
    #check power identity
    #print 'power1 is', (U_1_gpu.get()*Rho_gpu.get()+U_2_gpu.get()*Rho_gpu.get()).sum()*dtheta-(u_1_gpu.get()*rho_1_gpu.get()+u_2_gpu.get()*rho_2_gpu.get()).sum()*dx*dx
    #print 'power2 is', (U_1_gpu.get()*F_1_gpu.get()+U_2_gpu.get()*F_2_gpu.get()).sum()*dtheta-(u_1_gpu.get()*u_1_tilde_gpu.get()+u_2_gpu.get()*u_2_tilde_gpu.get()).sum()*dx*dx
    #check momentum conservation
    #print 'mom1 is', (u_1_gpu.get()*(rho_1_gpu.get()+rho_0)).sum()*dx*dx
    #print 'mom2 is', (u_2_gpu.get()*(rho_2_gpu.get()+rho_0)).sum()*dx*dx
    #print 'mom1 is', u_1_gpu.get().sum()*dx*dx
    #print 'mom2 is', u_2_gpu.get().sum()*dx*dx
    if PR_INT:
        vort_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), u_1_gpu.data, u_2_gpu.data, dx, N, vort_gpu.data)

    add_gpu(f_1_gpu, u_1_tilde_gpu, 1, f_1_gpu) 
    add_gpu(f_2_gpu, u_2_tilde_gpu, 1, f_2_gpu) 

    vc_laplace_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), u_1_gpu.data, u_2_gpu.data, mu_1_gpu.data, mu_2_gpu.data, dx, N, u_1_tilde_gpu.data, u_2_tilde_gpu.data) 

    add_gpu(f_1_gpu, u_1_tilde_gpu, 1.0, f_1_gpu) 
    add_gpu(f_2_gpu, u_2_tilde_gpu, 1.0, f_2_gpu) 

    ss_conv_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), u_1_gpu.data, u_2_gpu.data, dx, N, u_1_tilde_gpu.data, u_2_tilde_gpu.data)
    #zero_gpu(u_1_tilde_gpu) #turn off convection
    #zero_gpu(u_2_tilde_gpu)

    add_gpu(f_1_gpu, u_1_tilde_gpu, -rho_hat, f_1_gpu)
    add_gpu(f_2_gpu, u_2_tilde_gpu, -rho_hat, f_2_gpu)

    dal_term_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), u_1_gpu.data, u_2_gpu.data, u_1_o_gpu.data, u_2_o_gpu.data, u_1_tilde_gpu.data, u_2_tilde_gpu.data, rho_1_gpu.data, rho_2_gpu.data, dt, N, mu_1_gpu.data, mu_2_gpu.data)

    copy_gpu(u_1_gpu, u_1_o_gpu)
    copy_gpu(u_2_gpu, u_2_o_gpu)

    copy_gpu(f_1_gpu, u_1_tilde_gpu)
    copy_gpu(f_2_gpu, u_2_tilde_gpu)

    add_gpu(f_1_gpu, mu_1_gpu, -1.0, f_1_gpu)
    add_gpu(f_2_gpu, mu_2_gpu, -1.0, f_2_gpu)
    ##########################################################solve for u_tilde

    #make sure starting velocity is div free 
    add_gpu(f_1_gpu, u_1_gpu, rho_hat/dt, rho_1_extra_gpu)
    add_gpu(f_2_gpu, u_2_gpu, rho_hat/dt, rho_2_extra_gpu)
    div_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), rho_1_extra_gpu.data, rho_2_extra_gpu.data, dx, N, p_gpu.data)

    #div_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), f_1_gpu.data, f_2_gpu.data, dx, N, p_gpu.data)
    imag_gpu(p_gpu, p_tf_gpu)
    plan.execute(p_tf_gpu.data)

    L_adj_2_kernel(queue, (1,1), (1,1), L_tf_gpu.data)
    calc_p_gpu(p_tf_gpu, L_tf_gpu, p_tf_gpu)
    p_filt_kernel(queue, (1,1), (1,1), p_tf_gpu.data)

    plan.execute(p_tf_gpu.data, inverse=True)
    real_gpu(p_tf_gpu, p_gpu)
    #print 'p is', p_gpu.get().max()

    grad_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), p_gpu.data, dx, N, u_1_gpu.data, u_2_gpu.data)

    add_gpu(f_1_gpu, u_1_gpu, -1, f_1_gpu)
    add_gpu(f_2_gpu, u_2_gpu, -1, f_2_gpu)

    imag_gpu(f_1_gpu, f_1_tf_gpu)
    imag_gpu(f_2_gpu, f_2_tf_gpu)

    plan.execute(f_1_tf_gpu.data)
    plan.execute(f_2_tf_gpu.data)

    L_adj_1_kernel(queue, (1,1), (1,1), L_tf_gpu.data)
    solve_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), u_1_tf_gpu.data, u_2_tf_gpu.data, L_tf_gpu.data, f_1_tf_gpu.data, f_2_tf_gpu.data, dt, rho_hat, mu_0-mu_hat, mu_hat, N, f_1_tf_gpu.data, f_2_tf_gpu.data)

    #nyq_filt_kernel(queue, (N,1), (WGPSIZE_1, 1), f_1_tf_gpu.data, N)
    #nyq_filt_kernel(queue, (N,1), (WGPSIZE_1, 1), f_2_tf_gpu.data, N)

    zero_gpu(f_1_gpu)
    zero_gpu(f_2_gpu)

    #calc_f_gpu(dt*t_step, dx, xx_1_stag_gpu, xx_2_stag_gpu, xx_1_c_gpu, xx_2_c_gpu, f_1_gpu, f_2_gpu) 
    #scale_gpu(f_1_gpu, 0.5, f_1_gpu)
    #scale_gpu(f_2_gpu, 0.5, f_2_gpu)

    add_gpu(f_1_gpu, u_1_tilde_gpu, 0.5, f_1_gpu)
    add_gpu(f_2_gpu, u_2_tilde_gpu, 0.5, f_2_gpu)

    plan.execute(f_1_tf_gpu.data, inverse=True)
    plan.execute(f_2_tf_gpu.data, inverse=True)
    
    real_gpu(f_1_tf_gpu, u_1_tilde_gpu)
    real_gpu(f_2_tf_gpu, u_2_tilde_gpu)
    #print 'u_1 is', u_1_tilde_gpu.get().max()
    #print 'u_2 is', u_2_tilde_gpu.get().max()

    #copy_gpu(u_1_tilde_gpu,u_1_gpu)
    #copy_gpu(u_2_tilde_gpu,u_2_gpu)
    #copy_gpu(X_1_tilde_gpu,X_1_gpu)
    #copy_gpu(X_2_tilde_gpu,X_2_gpu)

    #return
    ############################################################

    ############################################################solve for u_soln
    zero_gpu(F_1_gpu)
    zero_gpu(F_2_gpu)

    update_routine(queue, firstn1_gpu, nextn1_gpu, firstn2_gpu, nextn2_gpu, X_1_gpu, X_2_gpu, N, dx)
    #update_routine(queue, firstn1_gpu, nextn1_gpu, firstn2_gpu, nextn2_gpu, X_1_tilde_gpu, X_2_tilde_gpu, N, dx)

    zero_gpu(u_1_gpu)
    zero_gpu(u_2_gpu)
    zero_gpu(mu_1_gpu)
    zero_gpu(mu_2_gpu)
    zero_gpu(rho_1_extra_gpu)
    zero_gpu(rho_2_extra_gpu)

    calc_lagF_teth_kernel(queue, (N_teth, 1), (b_WGPSIZE_4, 1), X_1_gpu.data, X_2_gpu.data, X_1_teth_gpu.data, X_2_teth_gpu.data, N_pre_teth, k_teth, dx_teth, F_1_gpu.data, F_2_gpu.data)

    spread_routine(queue, firstn1_gpu, nextn1_gpu, firstn2_gpu, nextn2_gpu, F_1_gpu, F_2_gpu, X_1_gpu, X_2_gpu, N, dx, dtheta, u_1_gpu, u_2_gpu)
    #spread_routine(queue, firstn1_gpu, nextn1_gpu, firstn2_gpu, nextn2_gpu, F_1_gpu, F_2_gpu, X_1_tilde_gpu, X_2_tilde_gpu, N, dx, dtheta, u_1_gpu, u_2_gpu)

    add_gpu(u_2_gpu, region_gpu, body_f*np.cos(wave_omega*dt*t_step), u_2_gpu)
    #add_gpu(u_2_gpu, region_gpu, body_f*np.cos(0.71*dt*t_step), u_2_gpu)

    calc_coefs_kernel(queue, (N_tri, 1), (b_WGPSIZE_3,1), X_1_tilde_gpu.data, X_2_tilde_gpu.data, TV_1_gpu.data, TV_2_gpu.data, TV_3_gpu.data, Rho_gpu.data, Mu_gpu.data, dx, N, rho_1_extra_gpu.data, rho_2_extra_gpu.data, mu_1_gpu.data, mu_2_gpu.data) 
    #calc_coefs_kernel(queue, (1, 1), (1,1), X_1_tilde_gpu.data, X_2_tilde_gpu.data, TV_1_gpu.data, TV_2_gpu.data, TV_3_gpu.data, Rho_gpu.data, Mu_gpu.data, dx, N, rho_1_extra_gpu.data, rho_2_extra_gpu.data, mu_1_gpu.data, mu_2_gpu.data) 

    zero_gpu(en_sum_gpu)
    calc_rho_kernel(queue, (N_tri, 1), (b_WGPSIZE_3, 1), X_1_tilde_gpu.data, X_2_tilde_gpu.data, TV_1_gpu.data, TV_2_gpu.data, TV_3_gpu.data, Rho_gpu.data, en_sum_gpu.data)

    rho_av = (en_sum_gpu.get().sum()+dx*dx*N*N)/(dx*dx*N*N)
    rho_av = (rho_0*N*N+rho_2_extra_gpu.get().sum())/(N*N) 
    #print 'rho_av is', rho_av
    #rho_av = 0
    #rho_av = np.mean(myrho_0*np.exp(-2*pi*xx_2_stag/C_h))

    #add_gpu(u_2_gpu, rho_0-rho_av+rho_2_extra_gpu, -grav, u_2_gpu)

    #print u_1_gpu.get()
    #print 'av force is', np.mean(u_2_gpu.get())

    add_gpu(rho_1_gpu, rho_1_extra_gpu, 1.0, rho_1_gpu)
    add_gpu(rho_2_gpu, rho_2_extra_gpu, 1.0, rho_2_gpu)
    #print 'av force is', np.mean(u_2_gpu.get())

    add_scalar_gpu(rho_2_extra_gpu, rho_0-rho_av, rho_2_extra_gpu)
    add_gpu(u_2_gpu, rho_2_extra_gpu, -grav, u_2_gpu)

    #add_scalar_gpu(mu_1_gpu, -mu_bar, mu_1_gpu)
    #add_scalar_gpu(mu_2_gpu, -mu_bar, mu_2_gpu)
    add_scalar_gpu(mu_1_gpu, mu_0-mu_hat, mu_1_gpu)
    add_scalar_gpu(mu_2_gpu, mu_0-mu_hat, mu_2_gpu)
    #add_scalar_gpu(rho_1_gpu, -rho_bar, rho_1_gpu)
    #add_scalar_gpu(rho_2_gpu, -rho_bar, rho_2_gpu)
    add_scalar_gpu(rho_1_gpu, rho_0-rho_hat, rho_1_gpu)
    add_scalar_gpu(rho_2_gpu, rho_0-rho_hat, rho_2_gpu)
    if SAVE:
        copy_gpu(mu_1_gpu+mu_0, mu_1_save_gpu)
        copy_gpu(mu_2_gpu+mu_0, mu_2_save_gpu)
        copy_gpu(rho_1_gpu+rho_0, rho_1_save_gpu)
        copy_gpu(rho_2_gpu+rho_0, rho_2_save_gpu)

    scale_gpu(rho_1_gpu, 0.5, rho_1_gpu)
    scale_gpu(rho_2_gpu, 0.5, rho_2_gpu)

    add_gpu(f_1_gpu, u_1_gpu, 0.5, f_1_gpu)
    add_gpu(f_2_gpu, u_2_gpu, 0.5, f_2_gpu)

    vc_laplace_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), u_1_tilde_gpu.data, u_2_tilde_gpu.data, mu_1_gpu.data, mu_2_gpu.data, dx, N, u_1_gpu.data, u_2_gpu.data) 

    add_gpu(f_1_gpu, u_1_gpu, 0.5, f_1_gpu)
    add_gpu(f_2_gpu, u_2_gpu, 0.5, f_2_gpu)

    ss_conv_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), u_1_tilde_gpu.data, u_2_tilde_gpu.data, dx, N, u_1_gpu.data, u_2_gpu.data)
    #zero_gpu(u_1_gpu) #turn off convection
    #zero_gpu(u_2_gpu)

    add_gpu(f_1_gpu, u_1_gpu, -rho_hat*0.5, f_1_gpu)
    add_gpu(f_2_gpu, u_2_gpu, -rho_hat*0.5, f_2_gpu)

    ss_conv_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), u_1_o_gpu.data, u_2_o_gpu.data, dx, N, mu_1_gpu.data, mu_2_gpu.data)
    #zero_gpu(mu_1_gpu) #turn off convection
    #zero_gpu(mu_2_gpu)

    add_gpu(u_1_gpu, mu_1_gpu, 1.0, u_1_gpu)
    add_gpu(u_2_gpu, mu_2_gpu, 1.0, u_2_gpu)

    scale_gpu(u_1_gpu, 0.5, u_1_gpu)
    scale_gpu(u_2_gpu, 0.5, u_2_gpu)

    dal_term_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), u_1_tilde_gpu.data, u_2_tilde_gpu.data, u_1_o_gpu.data, u_2_o_gpu.data, u_1_gpu.data, u_2_gpu.data, rho_1_gpu.data, rho_2_gpu.data, dt, N, u_1_gpu.data, u_2_gpu.data)

    add_gpu(f_1_gpu, u_1_gpu, -1.0, f_1_gpu)
    add_gpu(f_2_gpu, u_2_gpu, -1.0, f_2_gpu)

    #make sure starting velocity is div free 
    add_gpu(f_1_gpu, u_1_o_gpu, rho_hat/dt, rho_1_extra_gpu)
    add_gpu(f_2_gpu, u_2_o_gpu, rho_hat/dt, rho_2_extra_gpu)
    div_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), rho_1_extra_gpu.data, rho_2_extra_gpu.data, dx, N, p_gpu.data)

    #div_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), f_1_gpu.data, f_2_gpu.data, dx, N, p_gpu.data)
    imag_gpu(p_gpu, p_tf_gpu)
    plan.execute(p_tf_gpu.data)

    L_adj_2_kernel(queue, (1,1), (1,1), L_tf_gpu.data)
    calc_p_gpu(p_tf_gpu, L_tf_gpu, p_tf_gpu)
    p_filt_kernel(queue, (1,1), (1,1), p_tf_gpu.data)
    #nyq_filt_kernel(queue, (N,1), (WGPSIZE_1, 1), p_tf_gpu.data, N)
    plan.execute(p_tf_gpu.data, inverse=True)
    real_gpu(p_tf_gpu, p_gpu)
    #print 'p is', p_gpu.get().max()

    grad_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), p_gpu.data, dx, N, u_1_gpu.data, u_2_gpu.data)

    add_gpu(f_1_gpu, u_1_gpu, -1, f_1_gpu)
    add_gpu(f_2_gpu, u_2_gpu, -1, f_2_gpu)

    imag_gpu(f_1_gpu, f_1_tf_gpu)
    imag_gpu(f_2_gpu, f_2_tf_gpu)

    plan.execute(f_1_tf_gpu.data)
    plan.execute(f_2_tf_gpu.data)

    L_adj_1_kernel(queue, (1,1), (1,1), L_tf_gpu.data)
    solve_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), u_1_tf_gpu.data, u_2_tf_gpu.data, L_tf_gpu.data, f_1_tf_gpu.data, f_2_tf_gpu.data, dt, rho_hat, mu_0/2.0, mu_hat/2.0, N, f_1_tf_gpu.data, f_2_tf_gpu.data)

    #nyq_filt_kernel(queue, (N,1), (WGPSIZE_1, 1), f_1_tf_gpu.data, N)
    #nyq_filt_kernel(queue, (N,1), (WGPSIZE_1, 1), f_2_tf_gpu.data, N)

    plan.execute(f_1_tf_gpu.data, inverse=True)
    plan.execute(f_2_tf_gpu.data, inverse=True)
    
    real_gpu(f_1_tf_gpu, u_1_gpu)
    real_gpu(f_2_tf_gpu, u_2_gpu)

    interp_kernel(queue, (b_q+N_teth+b_add, 1), (b_WGPSIZE_2, 1), u_1_gpu.data, u_2_gpu.data, X_1_tilde_gpu.data, X_2_tilde_gpu.data, N, dx, F_1_gpu.data, F_2_gpu.data)
    #interp_kernel(queue, (b_q, 1), (b_WGPSIZE_2, 1), u_1_gpu.data, u_2_gpu.data, X_1_tilde_gpu.data, X_2_tilde_gpu.data, N, dx, F_1_gpu.data, F_2_gpu.data)
    add_gpu(U_1_gpu, F_1_gpu, 1, U_1_gpu)
    add_gpu(U_2_gpu, F_2_gpu, 1, U_2_gpu)
    add_gpu(X_1_gpu, U_1_gpu, dt/2, X_1_gpu)
    add_gpu(X_2_gpu, U_2_gpu, dt/2, X_2_gpu)
    update_routine(queue, firstn1_gpu, nextn1_gpu, firstn2_gpu, nextn2_gpu, X_1_gpu, X_2_gpu, N, dx)

    if SAVE:
        zero_gpu(mu_1_gpu)
        zero_gpu(mu_2_gpu)
        zero_gpu(rho_1_gpu)
        zero_gpu(rho_2_gpu)
        calc_coefs_kernel(queue, (N_tri, 1), (b_WGPSIZE_3,1), X_1_gpu.data, X_2_gpu.data, TV_1_gpu.data, TV_2_gpu.data, TV_3_gpu.data, Rho_gpu.data, Mu_gpu.data, dx, N, rho_1_gpu.data, rho_2_gpu.data, mu_1_gpu.data, mu_2_gpu.data) 
        copy_gpu(rho_1_gpu, rho_1_save_gpu)
        copy_gpu(rho_2_gpu, rho_2_save_gpu)
        add_scalar_gpu(rho_1_save_gpu, rho_0, rho_1_save_gpu)
        add_scalar_gpu(rho_2_save_gpu, rho_0, rho_2_save_gpu)
        copy_gpu(mu_1_gpu, mu_1_save_gpu)
        copy_gpu(mu_2_gpu, mu_2_save_gpu)
        add_scalar_gpu(mu_1_save_gpu, mu_0, mu_1_save_gpu)
        add_scalar_gpu(mu_2_save_gpu, mu_0, mu_2_save_gpu)
        #zero_gpu(F_1_gpu)
        #zero_gpu(F_2_gpu)
        #zero_gpu(f_1_gpu)
        #zero_gpu(f_2_gpu)
        #spread_routine(queue, firstn1_gpu, nextn1_gpu, firstn2_gpu, nextn2_gpu, F_1_gpu, F_2_gpu, X_1_gpu, X_2_gpu, N, dx, dtheta, f_1_gpu, f_2_gpu)
    ##############################################################
    return

def err(u_1, u_2, u_1_o, u_2_o, X_1, X_2, X_1_o, X_2_o, p, p_o, f_1, f_2, f_1_o, f_2_o):
    global dx, dtheta
    NN_o = u_1_o.shape[0]
    errinf1 = 0.
    errinf2 = 0.
    err11 = 0.
    err12 = 0.
    err21 = 0.
    err22 = 0.

    #NN_r_o = NN_o/32
    #NN_s_o = NN_o/2
    #NN_r = NN_o/16
    #NN_s = NN_o
    #NN_r_o = NN_o/8
    #NN_s_o = 25*NN_o/16
    #NN_r = NN_o/4
    #NN_s = 25*NN_o/8
    NN_r_o = NN_o/8
    NN_s_o = 50*NN_o/16
    NN_r = NN_o/4
    NN_s = 50*NN_o/8
    #NN_r_o = 3*NN_o/8
    #NN_s_o = 75*NN_o/16
    #NN_r = 3*NN_o/4
    #NN_s = 75*NN_o/8
    #X_1_o = X_1_o.reshape(NN_s_o,NN_r_o)
    #X_2_o = X_2_o.reshape(NN_s_o,NN_r_o)
    #X_1 = X_1.reshape(NN_s,NN_r)
    #X_2 = X_2.reshape(NN_s,NN_r)
    X_1_o = X_1_o.reshape(NN_r_o,NN_s_o).T
    X_2_o = X_2_o.reshape(NN_r_o,NN_s_o).T
    X_1 = X_1.reshape(NN_r,NN_s).T
    X_2 = X_2.reshape(NN_r,NN_s).T
    errinfX = 0.
    err1X = 0.
    err2X = 0.

    errinfp = 0
    err1p = 0
    err2p = 0

    vis_err_p = np.zeros((NN_o,NN_o))
    vis_err_u_1 = np.zeros((NN_o,NN_o))
    vis_err_u_2 = np.zeros((NN_o,NN_o))
    vis_err_f_1 = np.zeros((NN_o,NN_o))
    vis_err_f_2 = np.zeros((NN_o,NN_o))
    dx_o = 1.0/NN_o
    x_1_c_o = np.arange(-0.5+dx_o/2, 0.5+dx_o/2, dx_o)
    x_2_c_o = np.arange(-0.5+dx_o/2, 0.5+dx_o/2, dx_o)
    xx_1_c_o, xx_2_c_o = np.meshgrid(x_1_c_o, x_2_c_o)
    mesh3 = np.vstack((xx_1_c_o[np.newaxis],xx_2_c_o[np.newaxis],np.zeros((NN_o,NN_o))[np.newaxis]))

    for i in range(0,NN_o):
        for j in range(0,NN_o):
            vis_err_f_1[j,i] = np.abs(f_1_o[j,i]-0.5*(f_1[2*j,2*i]+f_1[2*j+1,2*i])) 
            vis_err_f_2[j,i] = np.abs(f_2_o[j,i]-0.5*(f_2[2*j,2*i]+f_2[2*j+1,2*i])) 
            vis_err_u_1[j,i] = np.abs(u_1_o[j,i]-0.5*(u_1[2*j,2*i]+u_1[2*j+1,2*i])) 
            if (np.abs(u_1_o[j,i]-0.5*(u_1[2*j,2*i]+u_1[2*j+1,2*i]))>errinf1):
                errinf1 = np.abs(u_1_o[j,i]-0.5*(u_1[2*j,2*i]+u_1[2*j+1,2*i]))
            vis_err_u_2[j,i] = np.abs(u_2_o[j,i]-0.5*(u_2[2*j,2*i]+u_2[2*j,2*i+1])) 
            if (np.abs(u_2_o[j,i]-0.5*(u_2[2*j,2*i]+u_2[2*j,2*i+1]))>errinf2):
                errinf2 = np.abs(u_2_o[j,i]-0.5*(u_2[2*j,2*i]+u_2[2*j,2*i+1]))
            err21 += np.abs(u_1_o[j,i]-0.5*(u_1[2*j,2*i]+u_1[2*j+1,2*i]))**2
            err22 += np.abs(u_2_o[j,i]-0.5*(u_2[2*j,2*i]+u_2[2*j,2*i+1]))**2
            err11 += np.abs(u_1_o[j,i]-0.5*(u_1[2*j,2*i]+u_1[2*j+1,2*i]))
            err12 += np.abs(u_2_o[j,i]-0.5*(u_2[2*j,2*i]+u_2[2*j,2*i+1]))
            vis_err_p[j,i] = np.abs(p_o[j,i]-0.25*(p[2*j,2*i]+p[2*j+1,2*i]+p[2*j,2*i+1]+p[2*j+1,2*i+1]))
            if (np.abs(p_o[j,i]-0.25*(p[2*j,2*i]+p[2*j+1,2*i]+p[2*j,2*i+1]+p[2*j+1,2*i+1]))>errinfp):
                errinfp = np.abs(p_o[j,i]-0.25*(p[2*j,2*i]+p[2*j+1,2*i]+p[2*j,2*i+1]+p[2*j+1,2*i+1]))
            err2p += np.abs(p_o[j,i]-0.25*(p[2*j,2*i]+p[2*j+1,2*i]+p[2*j,2*i+1]+p[2*j+1,2*i+1]))**2
            err1p += np.abs(p_o[j,i]-0.25*(p[2*j,2*i]+p[2*j+1,2*i]+p[2*j,2*i+1]+p[2*j+1,2*i+1])) 

    from os.path import join as pjoin
    vis_err_p = vis_err_p.reshape(NN_o,NN_o,1)
    file_name = "p-err-%04d.vts"%ii
    path_to_file = pjoin("vis_files", file_name)
    write_structured_grid(path_to_file, mesh3.reshape(3, NN_o, NN_o, 1),
    point_data=[("p-err", vis_err_p)])

    vis_err_f_1 = vis_err_f_1.reshape(NN_o,NN_o,1)
    file_name = "f-err_1-%04d.vts"%ii
    path_to_file = pjoin("vis_files", file_name) 
    write_structured_grid(path_to_file, mesh3.reshape(3, NN_o, NN_o, 1), 
    point_data=[("f-err_1", vis_err_f_1)])

    vis_err_f_2 = vis_err_f_2.reshape(NN_o,NN_o,1)
    file_name = "f-err_2-%04d.vts"%ii
    path_to_file = pjoin("vis_files", file_name) 
    write_structured_grid(path_to_file, mesh3.reshape(3, NN_o, NN_o, 1), 
    point_data=[("f-err_2", vis_err_f_2)])

    vis_err_u_1 = vis_err_u_1.reshape(NN_o,NN_o,1)
    file_name = "u-err_1-%04d.vts"%ii
    path_to_file = pjoin("vis_files", file_name) 
    write_structured_grid(path_to_file, mesh3.reshape(3, NN_o, NN_o, 1), 
    point_data=[("u-err_1", vis_err_u_1)])

    vis_err_u_2 = vis_err_u_2.reshape(NN_o,NN_o,1)
    file_name = "u-err_2-%04d.vts"%ii
    path_to_file = pjoin("vis_files", file_name) 
    write_structured_grid(path_to_file, mesh3.reshape(3, NN_o, NN_o, 1), 
    point_data=[("u-err_2", vis_err_u_2)])

    err21 = (dx/2)*np.sqrt(err21)
    err22 = (dx/2)*np.sqrt(err22)
    err11 = (dx**2/4)*err11
    err12 = (dx**2/4)*err12
    err2p = (dx/2)*np.sqrt(err2p)
    err1p = (dx**2/4)*err1p

    for i in range(0,NN_r_o):
        for j in range(0,NN_s_o):
            if (np.sqrt((X_1_o[j,i]-0.25*(X_1[2*j,2*i]+X_1[2*j+1,2*i]+X_1[2*j,2*i+1]+X_1[2*j+1,2*i+1]))**2
                +(X_2_o[j,i]-0.25*(X_2[2*j,2*i]+X_2[2*j+1,2*i]+X_2[2*j,2*i+1]+X_2[2*j+1,2*i+1]))**2)>errinfX):
                errinfX = np.sqrt((X_1_o[j,i]-0.25*(X_1[2*j,2*i]+X_1[2*j+1,2*i]+X_1[2*j,2*i+1]+X_1[2*j+1,2*i+1]))**2
                +(X_2_o[j,i]-0.25*(X_2[2*j,2*i]+X_2[2*j+1,2*i]+X_2[2*j,2*i+1]+X_2[2*j+1,2*i+1]))**2)
            err2X += np.sqrt((X_1_o[j,i]-0.25*(X_1[2*j,2*i]+X_1[2*j+1,2*i]+X_1[2*j,2*i+1]+X_1[2*j+1,2*i+1]))**2
                +(X_2_o[j,i]-0.25*(X_2[2*j,2*i]+X_2[2*j+1,2*i]+X_2[2*j,2*i+1]+X_2[2*j+1,2*i+1]))**2)**2
            err1X += np.sqrt((X_1_o[j,i]-0.25*(X_1[2*j,2*i]+X_1[2*j+1,2*i]+X_1[2*j,2*i+1]+X_1[2*j+1,2*i+1]))**2
                +(X_2_o[j,i]-0.25*(X_2[2*j,2*i]+X_2[2*j+1,2*i]+X_2[2*j,2*i+1]+X_2[2*j+1,2*i+1]))**2)

    err2X = np.sqrt(dtheta/4)*np.sqrt(err2X)
    err1X = (dtheta/4)*err1X

    return errinf1, errinf2, err11, err12, err21, err22, errinfX, err1X, err2X, errinfp, err1p, err2p 

ctx = cl.create_some_context()
queue = cl.CommandQueue(ctx)
Nyq_gpu = ElementwiseKernel(ctx,
	"double2 *u, double *filt, double2 *f",
	"f[i] = filt[i]*u[i]",
	"Nyq_gpu",
	preamble="""
        #pragma OPENCL EXTENSION cl_khr_fp64 : enable
	""")

calc_lagF_gpu = ElementwiseKernel(ctx,
	"double *X_1, double *X_2, "
	"double *rr_lag, unsigned *ip_s_lag, unsigned *im_s_lag, " 
	"double ds, "
	"double *F_1, double *F_2",
	"F_1[i] = 1./(ds*ds)*(1.+sin(2*pi*rr_lag[i]-pi/2))"
	   "*(X_1[ip_s_lag[i]]-2.*X_1[i]+X_1[im_s_lag[i]]); "
	"F_2[i] = 1./(ds*ds)*(1.+sin(2*pi*rr_lag[i]-pi/2))"
	   "*(X_2[ip_s_lag[i]]-2.*X_2[i]+X_2[im_s_lag[i]])",
	#"F_1[i] = 1./(ds*ds)"
	#   "*(X_1[ip_s_lag[i]]-2.*X_1[i]+X_1[im_s_lag[i]]); "
	#"F_2[i] = 1./(ds*ds)"
	#   "*(X_2[ip_s_lag[i]]-2.*X_2[i]+X_2[im_s_lag[i]])",
        #"F_1[i] = 0.;"
        #"F_2[i] = 0.",
	"calc_lagF_gpu",
	preamble="""
	#define pi (double)3.141592653589793
        #pragma OPENCL EXTENSION cl_khr_fp64 : enable
	""")

complex_mult_gpu = ElementwiseKernel(ctx,
	"double2 *u_1, double2 *u_2, double2 *f",
	"f[i] = complex_mul(u_1[i], u_2[i])",
	"complex_mult_gpu",
	preamble="""
	#define complex_ctr(x, y) (double2)(x, y)
	#define complex_mul(a, b) complex_ctr(-(a).y*(b).y+(a).x*(b).x, (a).y*(b).x+(a).x*(b).y)
        #pragma OPENCL EXTENSION cl_khr_fp64 : enable
	""")

mult_gpu = ElementwiseKernel(ctx,
	"double *u_1, double *u_2, double *f",
	"f[i] = u_1[i]*u_2[i]",
	"mult_gpu",
	preamble="""
        #pragma OPENCL EXTENSION cl_khr_fp64 : enable
	""")

add_gpu = ElementwiseKernel(ctx,
	"double *u_1, double *u_2, double a, double *f",
	"f[i] = u_1[i]+a*u_2[i]",
	"add_gpu",
	preamble="""
        #pragma OPENCL EXTENSION cl_khr_fp64 : enable
	""")
    
add_scalar_gpu = ElementwiseKernel(ctx,
	"double *u, double a, double *f",
	"f[i] = u[i]+a",
	"add_scalar_gpu",
	preamble="""
        #pragma OPENCL EXTENSION cl_khr_fp64 : enable
	""")
    
scale_gpu = ElementwiseKernel(ctx,
	"double *u, double a, double *f",
	"f[i] = a*u[i]",
	"scale_gpu",
	preamble="""
        #pragma OPENCL EXTENSION cl_khr_fp64 : enable
	""")

copy_gpu = ElementwiseKernel(ctx,
	"double *u, double *f",
	"f[i] = u[i]",
	"copy_gpu",
	preamble="""
        #pragma OPENCL EXTENSION cl_khr_fp64 : enable
	""")

copy_int_gpu = ElementwiseKernel(ctx,
	"unsigned *u, unsigned *f",
	"f[i] = u[i]",
	"copy_int_gpu",
	preamble="""
        #pragma OPENCL EXTENSION cl_khr_fp64 : enable
	""")

real_gpu = ElementwiseKernel(ctx,
	"double2 *x, double *z",
	"z[i] = x[i].x",
	"real_gpu",
	preamble="""
        #pragma OPENCL EXTENSION cl_khr_fp64 : enable
	""")

imag_gpu = ElementwiseKernel(ctx,
	"double *x, double2 *z",
	"z[i].x = x[i];"
	"z[i].y = 0",
	"imag_gpu",
	preamble="""
        #pragma OPENCL EXTENSION cl_khr_fp64 : enable
	""")

calc_f_gpu = ElementwiseKernel(ctx,
        "double t, double dx,"
        "double *xx_1_stag, double *xx_2_stag,"
        "double *xx_1_c, double *xx_2_c,"
        "double *f_1, double *f_2",
        "f_1[i] = 0.5*(2*cos(2*pi*xx_1_stag[i])*sin(4*pi*xx_2_c[i])*cos(t)"
           "-4*pi*sin(t)*sin(t)*sin(4*pi*xx_1_stag[i])"
           "+2*pi*cos(2*pi*xx_1_stag[i])*sin(2*pi*xx_2_c[i])*cos(t)"
           "+0.01*(1.0+0.5)*40*pi*pi*cos(2*pi*xx_1_stag[i])*sin(4*pi*xx_2_c[i])*sin(t)"
           "+16*pi*pi*0.5*0.01*cos(4*pi*xx_1_stag[i])*cos(2*pi*xx_2_c[i])*sin(4*pi*xx_2_c[i])*sin(t)"
           "+0.5*0.01*12*pi*pi*cos(2*pi*xx_1_stag[i])*cos(2*pi*xx_1_stag[i])"
           "*sin(t)*(sin(2*pi*xx_2_c[i])*cos(4*pi*xx_2_c[i])"
           "+2*cos(2*pi*xx_2_c[i])*sin(4*pi*xx_2_c[i])));"
        "f_2[i] = 0.5*(-sin(2*pi*xx_1_c[i])*cos(4*pi*xx_2_stag[i])*cos(t)"
           "-2*pi*sin(t)*sin(t)*sin(8*pi*xx_2_stag[i])"
           "+2*pi*sin(2*pi*xx_1_c[i])*cos(2*pi*xx_2_stag[i])*cos(t)"
           "-0.01*(1.0+0.5)*20*pi*pi*sin(2*pi*xx_1_c[i])*cos(4*pi*xx_2_stag[i])*sin(t)"
           "+12*pi*pi*0.5*0.01*sin(4*pi*xx_1_c[i])*cos(2*pi*xx_2_stag[i])*cos(4*pi*xx_2_stag[i])*sin(t)"
           "-8*pi*pi*0.5*0.01*sin(4*pi*xx_1_c[i])*sin(t)*(-sin(2*pi*xx_2_stag[i])*sin(4*pi*xx_2_stag[i])"
           "+2*cos(2*pi*xx_2_stag[i])*cos(4*pi*xx_2_stag[i])))",
        "calc_f_gpu",
        preamble="""
	#define pi (double)3.141592653589793
        #pragma OPENCL EXTENSION cl_khr_fp64 : enable
        """)

zero_gpu = ElementwiseKernel(ctx,
	"double *x",
	"x[i] = 0",
	"zero_gpu",
	preamble="""
        #pragma OPENCL EXTENSION cl_khr_fp64 : enable
	""")
 
null_imag_gpu = ElementwiseKernel(ctx,
	"double2 *x, double2 *z",
	"z[i].x = x[i].x;"
	"z[i].y = 0",
	"null_imag_gpu",
	preamble="""
        #pragma OPENCL EXTENSION cl_khr_fp64 : enable
	""")

calc_p_gpu = ElementwiseKernel(ctx,
	"double2 *w_tf, double2 *L_tf, double2 *p_tf",
	"p_tf[i] = complex_div(w_tf[i], L_tf[i])",
	"calc_p_gpu",
	preamble="""
	#define complex_ctr(x, y) (double2)(x, y)
	#define complex_div(a, b) 1/((b).x*(b).x+(b).y*(b).y)*complex_ctr((a).y*(b).y+(a).x*(b).x,(a).y*(b).x-(a).x*(b).y)
        #pragma OPENCL EXTENSION cl_khr_fp64 : enable
	""")

if False:
    solve_gpu = ElementwiseKernel(ctx,
    "double2 *u_1_tf, double2 *u_2_tf, double2 *L_tf,"
    "double2 *f_1_tf, double2 *f_2_tf, double2 *dp_1_tf,"
    "double2 *dp_2_tf, double dt, double rho, double mu,"
    "double2 *u_1_tf_tilde, double2 *u_2_tf_tilde",
    "double2 rhs_1 = complex_mul(1+mu*dt/(2*rho)*L_tf[i], u_1_tf[i])+dt/rho*(f_1_tf[i]-dp_1_tf[i]);"
    "double2 rhs_2 = complex_mul(1+mu*dt/(2*rho)*L_tf[i], u_2_tf[i])+dt/rho*(f_2_tf[i]-dp_2_tf[i]);"
    "double2 lhs = 1-mu*dt/(2*rho)*L_tf[i];"
    "u_1_tf_tilde[i] = complex_div(rhs_1, lhs);"
    "u_2_tf_tilde[i] = complex_div(rhs_2, lhs);",
    "solve_gpu",
    preamble="""
    #define complex_ctr(x, y) (double2)(x, y)
    #define complex_mul(a, b) complex_ctr(mad(-(a).y, (b).y, (a).x * (b).x), mad((a).y, (b).x, (a).x * (b).y))
    #define complex_div(a, b) 1/((b).x*(b).x+(b).y*(b).y)*complex_ctr((a).y*(b).y+(a).x*(b).x,(a).y*(b).x-(a).x*(b).y)
    """)

ib_prg = cl.Program(ctx, """
    //#pragma OPENCL EXTENSION cl_amd_printf : enable
    #pragma OPENCL EXTENSION cl_khr_fp64 : enable
    //#pragma OPENCL EXTENSION cl_amd_fp64 : enable
    __kernel void interp(__global const double *u_1, __global const double *u_2, __global const double *X_1, __global const double *X_2, const unsigned N, const double dx, __global double *U_1, __global double *U_2)
    {
      unsigned k = get_global_id(0);
      double s_1_stag = (X_1[k]+0.5)/dx;
      double s_2_stag = (X_2[k]+0.5)/dx;       
      int i_1_stag = convert_int(floor(s_1_stag));
      int i_2_stag = convert_int(floor(s_2_stag));
      double r_1_stag = s_1_stag-i_1_stag;
      double r_2_stag = s_2_stag-i_2_stag;
      double s_1_c = (X_1[k]+0.5-dx/2)/dx;
      double s_2_c = (X_2[k]+0.5-dx/2)/dx;        
      int i_1_c = convert_int(floor(s_1_c));
      int i_2_c = convert_int(floor(s_2_c));
      double r_1_c = s_1_c-i_1_c;
      double r_2_c = s_2_c-i_2_c;
      
      int ind_1x, ind_1y, ind_1z, ind_1w;
      int ind_2x, ind_2y, ind_2z, ind_2w;
      ind_1x = (i_1_stag-1 & N-1);
      ind_1y = (i_1_stag & N-1);
      ind_1z = (i_1_stag+1 & N-1);
      ind_1w = (i_1_stag+2 & N-1);

      ind_2x = (i_2_c-1 & N-1);
      ind_2y = (i_2_c & N-1);
      ind_2z = (i_2_c+1 & N-1);
      ind_2w = (i_2_c+2 & N-1);

      double sb_tests0, sb_tests1, sb_tests2, sb_tests3;
      double sb_tests4, sb_tests5, sb_tests6, sb_tests7;
      double sb_tests8, sb_tests9, sb_testsa, sb_testsb;
      double sb_testsc, sb_testsd, sb_testse, sb_testsf;

      sb_tests0 = u_1[ind_1x+N*ind_2x];
      sb_tests1 = u_1[ind_1y+N*ind_2x];
      sb_tests2 = u_1[ind_1z+N*ind_2x];
      sb_tests3 = u_1[ind_1w+N*ind_2x];
      sb_tests4 = u_1[ind_1x+N*ind_2y];
      sb_tests5 = u_1[ind_1y+N*ind_2y];
      sb_tests6 = u_1[ind_1z+N*ind_2y];
      sb_tests7 = u_1[ind_1w+N*ind_2y];
      sb_tests8 = u_1[ind_1x+N*ind_2z];
      sb_tests9 = u_1[ind_1y+N*ind_2z];
      sb_testsa = u_1[ind_1z+N*ind_2z];
      sb_testsb = u_1[ind_1w+N*ind_2z];
      sb_testsc = u_1[ind_1x+N*ind_2w];
      sb_testsd = u_1[ind_1y+N*ind_2w];
      sb_testse = u_1[ind_1z+N*ind_2w];
      sb_testsf = u_1[ind_1w+N*ind_2w];

      double q1 = sqrt(1.+4.0*r_1_stag*(1-r_1_stag)); 
      double q2 = sqrt(1.+4.0*r_2_c*(1-r_2_c));
      double v1x, v1y, v1z, v1w;
      double v2x, v2y, v2z, v2w;

      v1x = 3.-2.*r_1_stag-q1; 
      v1y = 3.-2.*r_1_stag+q1; 
      v1z = 1.+2.*r_1_stag+q1; 
      v1w = 1.+2.*r_1_stag-q1; 
      v1x = v1x/8.;
      v1y = v1y/8.;
      v1z = v1z/8.;
      v1w = v1w/8.;

      v2x = 3.-2.*r_2_c-q2; 
      v2y = 3.-2.*r_2_c+q2; 
      v2z = 1.+2.*r_2_c+q2; 
      v2w = 1.+2.*r_2_c-q2; 
      v2x = v2x/8.;
      v2y = v2y/8.;
      v2z = v2z/8.;
      v2w = v2w/8.;
      
      double ws0, ws1, ws2, ws3;
      double ws4, ws5, ws6, ws7;
      double ws8, ws9, wsa, wsb;
      double wsc, wsd, wse, wsf;
      ws0 = v1x*v2x;
      ws1 = v1y*v2x;
      ws2 = v1z*v2x;
      ws3 = v1w*v2x;
      ws4 = v1x*v2y;
      ws5 = v1y*v2y;
      ws6 = v1z*v2y;
      ws7 = v1w*v2y;
      ws8 = v1x*v2z;
      ws9 = v1y*v2z;
      wsa = v1z*v2z;
      wsb = v1w*v2z;
      wsc = v1x*v2w;
      wsd = v1y*v2w;
      wse = v1z*v2w;
      wsf = v1w*v2w;
      
      ws0 = ws0*sb_tests0;
      ws1 = ws1*sb_tests1;
      ws2 = ws2*sb_tests2;
      ws3 = ws3*sb_tests3;
      ws4 = ws4*sb_tests4;
      ws5 = ws5*sb_tests5;
      ws6 = ws6*sb_tests6;
      ws7 = ws7*sb_tests7;
      ws8 = ws8*sb_tests8;
      ws9 = ws9*sb_tests9;
      wsa = wsa*sb_testsa;
      wsb = wsb*sb_testsb;
      wsc = wsc*sb_testsc;
      wsd = wsd*sb_testsd;
      wse = wse*sb_testse;
      wsf = wsf*sb_testsf;

      U_1[k] = ws0+ws1+ws2+ws3
	       +ws4+ws5+ws6+ws7
	       +ws8+ws9+wsa+wsb
	       +wsc+wsd+wse+wsf;

      ind_1x = (i_1_c-1 & N-1);
      ind_1y = (i_1_c & N-1);
      ind_1z = (i_1_c+1 & N-1);
      ind_1w = (i_1_c+2 & N-1);

      ind_2x = (i_2_stag-1 & N-1);
      ind_2y = (i_2_stag & N-1);
      ind_2z = (i_2_stag+1 & N-1);
      ind_2w = (i_2_stag+2 & N-1);

      sb_tests0 = u_2[ind_1x+N*ind_2x];
      sb_tests1 = u_2[ind_1y+N*ind_2x];
      sb_tests2 = u_2[ind_1z+N*ind_2x];
      sb_tests3 = u_2[ind_1w+N*ind_2x];
      sb_tests4 = u_2[ind_1x+N*ind_2y];
      sb_tests5 = u_2[ind_1y+N*ind_2y];
      sb_tests6 = u_2[ind_1z+N*ind_2y];
      sb_tests7 = u_2[ind_1w+N*ind_2y];
      sb_tests8 = u_2[ind_1x+N*ind_2z];
      sb_tests9 = u_2[ind_1y+N*ind_2z];
      sb_testsa = u_2[ind_1z+N*ind_2z];
      sb_testsb = u_2[ind_1w+N*ind_2z];
      sb_testsc = u_2[ind_1x+N*ind_2w];
      sb_testsd = u_2[ind_1y+N*ind_2w];
      sb_testse = u_2[ind_1z+N*ind_2w];
      sb_testsf = u_2[ind_1w+N*ind_2w];

      q1 = sqrt(1.+4.0*r_1_c*(1-r_1_c)); 
      q2 = sqrt(1.+4.0*r_2_stag*(1-r_2_stag));

      v1x = 3.-2.*r_1_c-q1; 
      v1y = 3.-2.*r_1_c+q1; 
      v1z = 1.+2.*r_1_c+q1; 
      v1w = 1.+2.*r_1_c-q1; 
      v1x = v1x/8.;
      v1y = v1y/8.;
      v1z = v1z/8.;
      v1w = v1w/8.;

      v2x = 3.-2.*r_2_stag-q2; 
      v2y = 3.-2.*r_2_stag+q2; 
      v2z = 1.+2.*r_2_stag+q2; 
      v2w = 1.+2.*r_2_stag-q2; 
      v2x = v2x/8.;
      v2y = v2y/8.;
      v2z = v2z/8.;
      v2w = v2w/8.;
      
      ws0 = v1x*v2x;
      ws1 = v1y*v2x;
      ws2 = v1z*v2x;
      ws3 = v1w*v2x;
      ws4 = v1x*v2y;
      ws5 = v1y*v2y;
      ws6 = v1z*v2y;
      ws7 = v1w*v2y;
      ws8 = v1x*v2z;
      ws9 = v1y*v2z;
      wsa = v1z*v2z;
      wsb = v1w*v2z;
      wsc = v1x*v2w;
      wsd = v1y*v2w;
      wse = v1z*v2w;
      wsf = v1w*v2w;
      
      ws0 = ws0*sb_tests0;
      ws1 = ws1*sb_tests1;
      ws2 = ws2*sb_tests2;
      ws3 = ws3*sb_tests3;
      ws4 = ws4*sb_tests4;
      ws5 = ws5*sb_tests5;
      ws6 = ws6*sb_tests6;
      ws7 = ws7*sb_tests7;
      ws8 = ws8*sb_tests8;
      ws9 = ws9*sb_tests9;
      wsa = wsa*sb_testsa;
      wsb = wsb*sb_testsb;
      wsc = wsc*sb_testsc;
      wsd = wsd*sb_testsd;
      wse = wse*sb_testse;
      wsf = wsf*sb_testsf;

      U_2[k] = ws0+ws1+ws2+ws3
	       +ws4+ws5+ws6+ws7
	       +ws8+ws9+wsa+wsb
	       +wsc+wsd+wse+wsf;
    }
    __kernel void spread(__global const int *firstn1, __global const int *nextn1, __global const int *firstn2, __global const int *nextn2, __global const double *F_1, __global const double *F_2, __global const double *X_1, __global const double *X_2, const unsigned N, const double dx, const double dtheta, const unsigned shift0, const unsigned shift1, __global double *f_1, __global double *f_2)
    {
      unsigned i = 4*get_global_id(0)+shift0;
      unsigned j = 4*get_global_id(1)+shift1;

      int k1 = firstn1[i+N*j]; 
      int k2 = firstn2[i+N*j]; 
      if (k1 == -1 && k2 == -1){
      return;
      }

      double c = dtheta/(dx*dx);
      int ind_1x, ind_1y, ind_1z, ind_1w;
      int ind_2x, ind_2y, ind_2z, ind_2w;

      double f_10, f_11, f_12, f_13;
      double f_14, f_15, f_16, f_17;
      double f_18, f_19, f_1a, f_1b;
      double f_1c, f_1d, f_1e, f_1f;

      double f_20, f_21, f_22, f_23;
      double f_24, f_25, f_26, f_27;
      double f_28, f_29, f_2a, f_2b;
      double f_2c, f_2d, f_2e, f_2f;

      double v1x, v1y, v1z, v1w;
      double v2x, v2y, v2z, v2w;

      double q1, q2;

      double ws0, ws1, ws2, ws3;
      double ws4, ws5, ws6, ws7;
      double ws8, ws9, wsa, wsb;
      double wsc, wsd, wse, wsf;

      if (k1 !=-1){
      ind_1x = (i-1 & N-1);
      ind_1y = (i & N-1);
      ind_1z = (i+1 & N-1);
      ind_1w = (i+2 & N-1);

      ind_2x = (j-1 & N-1);
      ind_2y = (j & N-1);
      ind_2z = (j+1 & N-1);
      ind_2w = (j+2 & N-1);

      f_10=f_1[ind_1x+N*ind_2x];
      f_11=f_1[ind_1y+N*ind_2x];
      f_12=f_1[ind_1z+N*ind_2x];
      f_13=f_1[ind_1w+N*ind_2x];
      f_14=f_1[ind_1x+N*ind_2y];
      f_15=f_1[ind_1y+N*ind_2y];
      f_16=f_1[ind_1z+N*ind_2y];
      f_17=f_1[ind_1w+N*ind_2y];
      f_18=f_1[ind_1x+N*ind_2z];
      f_19=f_1[ind_1y+N*ind_2z];
      f_1a=f_1[ind_1z+N*ind_2z];
      f_1b=f_1[ind_1w+N*ind_2z];
      f_1c=f_1[ind_1x+N*ind_2w];
      f_1d=f_1[ind_1y+N*ind_2w];
      f_1e=f_1[ind_1z+N*ind_2w];
      f_1f=f_1[ind_1w+N*ind_2w];
      }

      while(k1 != -1)
      {
      double s_1_stag = (X_1[k1]+0.5)/dx;
      int i_1_stag = convert_int(floor(s_1_stag));
      double r_1_stag = s_1_stag-i_1_stag;
      double s_2_c = (X_2[k1]+0.5-dx/2)/dx;        
      int i_2_c = convert_int(floor(s_2_c));
      double r_2_c = s_2_c-i_2_c;
      
      q1 = sqrt(1.+4.0*r_1_stag*(1-r_1_stag)); 
      q2 = sqrt(1.+4.0*r_2_c*(1-r_2_c));

      v1x = 3.-2.*r_1_stag-q1; 
      v1y = 3.-2.*r_1_stag+q1; 
      v1z = 1.+2.*r_1_stag+q1; 
      v1w = 1.+2.*r_1_stag-q1; 
      v1x = v1x/8.;
      v1y = v1y/8.;
      v1z = v1z/8.;
      v1w = v1w/8.;

      v2x = 3.-2.*r_2_c-q2; 
      v2y = 3.-2.*r_2_c+q2; 
      v2z = 1.+2.*r_2_c+q2; 
      v2w = 1.+2.*r_2_c-q2; 
      v2x = v2x/8.;
      v2y = v2y/8.;
      v2z = v2z/8.;
      v2w = v2w/8.;
      
      ws0 = v1x*v2x;
      ws1 = v1y*v2x;
      ws2 = v1z*v2x;
      ws3 = v1w*v2x;
      ws4 = v1x*v2y;
      ws5 = v1y*v2y;
      ws6 = v1z*v2y;
      ws7 = v1w*v2y;
      ws8 = v1x*v2z;
      ws9 = v1y*v2z;
      wsa = v1z*v2z;
      wsb = v1w*v2z;
      wsc = v1x*v2w;
      wsd = v1y*v2w;
      wse = v1z*v2w;
      wsf = v1w*v2w;
      
      double my_F_1 = F_1[k1];
      f_10 += c*my_F_1*ws0;
      f_11 += c*my_F_1*ws1;
      f_12 += c*my_F_1*ws2;
      f_13 += c*my_F_1*ws3;
      f_14 += c*my_F_1*ws4;
      f_15 += c*my_F_1*ws5;
      f_16 += c*my_F_1*ws6;
      f_17 += c*my_F_1*ws7;
      f_18 += c*my_F_1*ws8;
      f_19 += c*my_F_1*ws9;
      f_1a += c*my_F_1*wsa;
      f_1b += c*my_F_1*wsb;
      f_1c += c*my_F_1*wsc;
      f_1d += c*my_F_1*wsd;
      f_1e += c*my_F_1*wse;
      f_1f += c*my_F_1*wsf;

      k1 = nextn1[k1];
      if (k1 == -1){
      f_1[ind_1x+N*ind_2x]=f_10;
      f_1[ind_1y+N*ind_2x]=f_11;
      f_1[ind_1z+N*ind_2x]=f_12;
      f_1[ind_1w+N*ind_2x]=f_13;
      f_1[ind_1x+N*ind_2y]=f_14;
      f_1[ind_1y+N*ind_2y]=f_15;
      f_1[ind_1z+N*ind_2y]=f_16;
      f_1[ind_1w+N*ind_2y]=f_17;
      f_1[ind_1x+N*ind_2z]=f_18;
      f_1[ind_1y+N*ind_2z]=f_19;
      f_1[ind_1z+N*ind_2z]=f_1a;
      f_1[ind_1w+N*ind_2z]=f_1b;
      f_1[ind_1x+N*ind_2w]=f_1c;
      f_1[ind_1y+N*ind_2w]=f_1d;
      f_1[ind_1z+N*ind_2w]=f_1e;
      f_1[ind_1w+N*ind_2w]=f_1f;
      }
      }
      
      if (k2 !=-1){
      ind_1x = (i-1 & N-1);
      ind_1y = (i & N-1);
      ind_1z = (i+1 & N-1);
      ind_1w = (i+2 & N-1);

      ind_2x = (j-1 & N-1);
      ind_2y = (j & N-1);
      ind_2z = (j+1 & N-1);
      ind_2w = (j+2 & N-1);

      f_20=f_2[ind_1x+N*ind_2x];
      f_21=f_2[ind_1y+N*ind_2x];
      f_22=f_2[ind_1z+N*ind_2x];
      f_23=f_2[ind_1w+N*ind_2x];
      f_24=f_2[ind_1x+N*ind_2y];
      f_25=f_2[ind_1y+N*ind_2y];
      f_26=f_2[ind_1z+N*ind_2y];
      f_27=f_2[ind_1w+N*ind_2y];
      f_28=f_2[ind_1x+N*ind_2z];
      f_29=f_2[ind_1y+N*ind_2z];
      f_2a=f_2[ind_1z+N*ind_2z];
      f_2b=f_2[ind_1w+N*ind_2z];
      f_2c=f_2[ind_1x+N*ind_2w];
      f_2d=f_2[ind_1y+N*ind_2w];
      f_2e=f_2[ind_1z+N*ind_2w];
      f_2f=f_2[ind_1w+N*ind_2w];
      }
      
      while (k2 != -1){
      double s_1_c = (X_1[k2]+0.5-dx/2)/dx;
      int i_1_c = convert_int(floor(s_1_c));
      double r_1_c = s_1_c-i_1_c;
      double s_2_stag = (X_2[k2]+0.5)/dx;        
      int i_2_stag = convert_int(floor(s_2_stag));
      double r_2_stag = s_2_stag-i_2_stag;

      q1 = sqrt(1.+4.0*r_1_c*(1-r_1_c)); 
      q2 = sqrt(1.+4.0*r_2_stag*(1-r_2_stag));

      v1x = 3.-2.*r_1_c-q1; 
      v1y = 3.-2.*r_1_c+q1; 
      v1z = 1.+2.*r_1_c+q1; 
      v1w = 1.+2.*r_1_c-q1; 
      v1x = v1x/8.;
      v1y = v1y/8.;
      v1z = v1z/8.;
      v1w = v1w/8.;

      v2x = 3.-2.*r_2_stag-q2; 
      v2y = 3.-2.*r_2_stag+q2; 
      v2z = 1.+2.*r_2_stag+q2; 
      v2w = 1.+2.*r_2_stag-q2; 
      v2x = v2x/8.;
      v2y = v2y/8.;
      v2z = v2z/8.;
      v2w = v2w/8.;
      
      ws0 = v1x*v2x;
      ws1 = v1y*v2x;
      ws2 = v1z*v2x;
      ws3 = v1w*v2x;
      ws4 = v1x*v2y;
      ws5 = v1y*v2y;
      ws6 = v1z*v2y;
      ws7 = v1w*v2y;
      ws8 = v1x*v2z;
      ws9 = v1y*v2z;
      wsa = v1z*v2z;
      wsb = v1w*v2z;
      wsc = v1x*v2w;
      wsd = v1y*v2w;
      wse = v1z*v2w;
      wsf = v1w*v2w;
      
      double my_F_2 = F_2[k2];
      f_20 += c*my_F_2*ws0;
      f_21 += c*my_F_2*ws1;
      f_22 += c*my_F_2*ws2;
      f_23 += c*my_F_2*ws3;
      f_24 += c*my_F_2*ws4;
      f_25 += c*my_F_2*ws5;
      f_26 += c*my_F_2*ws6;
      f_27 += c*my_F_2*ws7;
      f_28 += c*my_F_2*ws8;
      f_29 += c*my_F_2*ws9;
      f_2a += c*my_F_2*wsa;
      f_2b += c*my_F_2*wsb;
      f_2c += c*my_F_2*wsc;
      f_2d += c*my_F_2*wsd;
      f_2e += c*my_F_2*wse;
      f_2f += c*my_F_2*wsf;

      k2 = nextn2[k2];
      if (k2 == -1){
      f_2[ind_1x+N*ind_2x]=f_20;
      f_2[ind_1y+N*ind_2x]=f_21;
      f_2[ind_1z+N*ind_2x]=f_22;
      f_2[ind_1w+N*ind_2x]=f_23;
      f_2[ind_1x+N*ind_2y]=f_24;
      f_2[ind_1y+N*ind_2y]=f_25;
      f_2[ind_1z+N*ind_2y]=f_26;
      f_2[ind_1w+N*ind_2y]=f_27;
      f_2[ind_1x+N*ind_2z]=f_28;
      f_2[ind_1y+N*ind_2z]=f_29;
      f_2[ind_1z+N*ind_2z]=f_2a;
      f_2[ind_1w+N*ind_2z]=f_2b;
      f_2[ind_1x+N*ind_2w]=f_2c;
      f_2[ind_1y+N*ind_2w]=f_2d;
      f_2[ind_1z+N*ind_2w]=f_2e;
      f_2[ind_1w+N*ind_2w]=f_2f;
      }
      }
    }
    __kernel void updaten(__global int *firstn1, __global int *nextn1, __global int *firstn2, __global int *nextn2, __global const double *X_1, __global const double *X_2, const unsigned N, const double dx, const unsigned shift0, const unsigned shift1)
    {
      unsigned i = 4*get_global_id(0)+shift0;
      unsigned j = 4*get_global_id(1)+shift1;

      int k1 = firstn1[i+N*j]; 
      int k2 = firstn2[i+N*j]; 
      if (k1 == -1 && k2 == -1){
      return;
      }

      uint ind_1, ind_2;
      int nextnold;

      int nprev1 = -1;
      while(k1 != -1)
      {
      double s_1_stag = (X_1[k1]+0.5)/dx;
      int i_1_stag = convert_int(floor(s_1_stag));
      double s_2_c = (X_2[k1]+0.5-dx/2)/dx;        
      int i_2_c = convert_int(floor(s_2_c));
       
      ind_1 = (i_1_stag & N-1);
      ind_2 = (i_2_c & N-1);

      if (ind_1 == i && ind_2 == j){
      nprev1 = k1;
      k1 = nextn1[k1];
      }
      else {
      nextnold = nextn1[k1];
      nextn1[k1] = firstn1[ind_1+N*ind_2];
      firstn1[ind_1+N*ind_2] = k1;
      k1 = nextnold;
          if (nprev1 == -1) {
          firstn1[i+N*j] = nextnold;
          }
          else {
          nextn1[nprev1] = nextnold;
          }
      }
      }

      int nprev2 = -1;
      while(k2 != -1)
      {
      double s_2_stag = (X_2[k2]+0.5)/dx;        
      int i_2_stag = convert_int(floor(s_2_stag));
      double s_1_c = (X_1[k2]+0.5-dx/2)/dx;
      int i_1_c = convert_int(floor(s_1_c));
      
      ind_1 = (i_1_c & N-1);
      ind_2 = (i_2_stag & N-1);

      if (ind_1 == i && ind_2 == j){
      nprev2 = k2;
      k2 = nextn2[k2];
      }
      else {
      nextnold = nextn2[k2];
      nextn2[k2] = firstn2[ind_1+N*ind_2];
      firstn2[ind_1+N*ind_2] = k2;
      k2 = nextnold;
          if (nprev2 == -1) {
          firstn2[i+N*j] = nextnold;
          }
          else {
          nextn2[nprev2] = nextnold;
          }
      }
      }
    }
    __kernel void calc_lagF_teth(__global const double *X_1, __global const double *X_2, __global const double *X_1_teth, __global const double *X_2_teth, const unsigned b_teth, const double k_teth, const double dr, __global double *F_1, __global double *F_2)
    { 
    unsigned i = get_global_id(0);

    F_1[i+b_teth] += -k_teth/dr*(X_1[i+b_teth]-X_1_teth[i]); 
    F_2[i+b_teth] += -k_teth/dr*(X_2[i+b_teth]-X_2_teth[i]); 
    }
    __kernel void calc_lagF(__global const double *X_1, __global const double *X_2, __global const unsigned *i_t, __global const unsigned *ip_t, __global const unsigned *im_t, const double k_s, const double k_b, const double ds, __global double *F_1, __global double *F_2) 
    {
    unsigned i = get_global_id(0);

    F_1[i] = 
             -k_s*((sqrt((X_1[i_t[i]]-X_1[im_t[i]])*(X_1[i_t[i]]-X_1[im_t[i]])+(X_2[i_t[i]]-X_2[im_t[i]])*(X_2[i_t[i]]-X_2[im_t[i]]))/ds-1)         
             *(X_1[i_t[i]]-X_1[im_t[i]])/sqrt((X_1[i_t[i]]-X_1[im_t[i]])*(X_1[i_t[i]]-X_1[im_t[i]])+(X_2[i_t[i]]-X_2[im_t[i]])*(X_2[i_t[i]]-X_2[im_t[i]]))
             -(sqrt((X_1[ip_t[i]]-X_1[i_t[i]])*(X_1[ip_t[i]]-X_1[i_t[i]])+(X_2[ip_t[i]]-X_2[i_t[i]])*(X_2[ip_t[i]]-X_2[i_t[i]]))/ds-1)
             *(X_1[ip_t[i]]-X_1[i_t[i]])/sqrt((X_1[ip_t[i]]-X_1[i_t[i]])*(X_1[ip_t[i]]-X_1[i_t[i]])+(X_2[ip_t[i]]-X_2[i_t[i]])*(X_2[ip_t[i]]-X_2[i_t[i]])))
             -k_b*(X_1[ip_t[ip_t[i]]]-4*X_1[ip_t[i]]+X_1[im_t[im_t[i]]]-4*X_1[im_t[i]]+6*X_1[i_t[i]])/(ds*ds*ds);
    F_2[i] = 
             -k_s*((sqrt((X_1[i_t[i]]-X_1[im_t[i]])*(X_1[i_t[i]]-X_1[im_t[i]])+(X_2[i_t[i]]-X_2[im_t[i]])*(X_2[i_t[i]]-X_2[im_t[i]]))/ds-1)         
             *(X_2[i_t[i]]-X_2[im_t[i]])/sqrt((X_1[i_t[i]]-X_1[im_t[i]])*(X_1[i_t[i]]-X_1[im_t[i]])+(X_2[i_t[i]]-X_2[im_t[i]])*(X_2[i_t[i]]-X_2[im_t[i]]))
             -(sqrt((X_1[ip_t[i]]-X_1[i_t[i]])*(X_1[ip_t[i]]-X_1[i_t[i]])+(X_2[ip_t[i]]-X_2[i_t[i]])*(X_2[ip_t[i]]-X_2[i_t[i]]))/ds-1)
             *(X_2[ip_t[i]]-X_2[i_t[i]])/sqrt((X_1[ip_t[i]]-X_1[i_t[i]])*(X_1[ip_t[i]]-X_1[i_t[i]])+(X_2[ip_t[i]]-X_2[i_t[i]])*(X_2[ip_t[i]]-X_2[i_t[i]])))
             -k_b*(X_2[ip_t[ip_t[i]]]-4*X_2[ip_t[i]]+X_2[im_t[im_t[i]]]-4*X_2[im_t[i]]+6*X_2[i_t[i]])/(ds*ds*ds);
    }
    __kernel void calc_coefs_retriang(__global const double *X_1, __global const double *X_2, __global const unsigned *TV_1, __global const unsigned *TV_2, __global const unsigned *TV_3, __global const double *Rho, __global const double *Mu, const unsigned N_inds, const unsigned N_pre, __global double *tmpRho, __global double *tmpMu)
    { 
    unsigned i = get_global_id(0);

    //possible improvement--bring Rho, Mu, TV into local memory

    double X_11, X_12;
    double X_21, X_22;
    double X_31, X_32;

    double mn1, mn2;
    double mx1, mx2;

    double gpt1, gpt2;

    double p121, p122;
    double p231, p232;
    double p311, p312;

    double chk121, chk122;
    double chk231, chk232;
    double chk311, chk312;

    double dot1, dot2, dot3;

    X_11 = X_1[TV_1[i]];
    X_12 = X_2[TV_1[i]];

    X_21 = X_1[TV_2[i]];
    X_22 = X_2[TV_2[i]];

    X_31 = X_1[TV_3[i]];
    X_32 = X_2[TV_3[i]];

    //find candidate grid points
    mn1 = fmin(X_11,X_21);
    mn1 = fmin(mn1,X_31);

    mn2 = fmin(X_12,X_22);
    mn2 = fmin(mn2,X_32);

    mx1 = fmax(X_11,X_21);
    mx1 = fmax(mx1,X_31);

    mx2 = fmax(X_12,X_22);
    mx2 = fmax(mx2,X_32);

    for (int k = 0; k < N_inds; k++) {
    //check whether candidate points lie inside triangle

       gpt1 = X_1[N_pre+k];
       gpt2 = X_2[N_pre+k];

       p121 = gpt1 - mn1 - floor(gpt1-mn1);
       p122 = gpt2 - mn2 - floor(gpt2-mn2);

       if (p121 > mx1-mn1 | p122 > mx2-mn2) {
           continue;
       }
       else {
           gpt1 = gpt1 - floor(gpt1-mn1);    
           gpt2 = gpt2 - floor(gpt2-mn2);    
       }

       //find perp vectors to half-spaces
       p121 = -(X_22-X_12);
       p122 = X_21-X_11;
       p231 = -(X_32-X_22);
       p232 = X_31-X_21;
       p311 = -(X_12-X_32);
       p312 = X_11-X_31;

       //find vectors to dot against
       chk121 = gpt1 - X_11;
       chk122 = gpt2 - X_12;
       chk231 = gpt1 - X_21;
       chk232 = gpt2 - X_22;
       chk311 = gpt1 - X_31;
       chk312 = gpt2 - X_32;

       dot1 = p121*chk121+p122*chk122;
       dot2 = p231*chk231+p232*chk232;
       dot3 = p311*chk311+p312*chk312;

       //constant interpolation
       //if (dot1 >= 0 && dot2 >= 0 && dot3 >= 0) {
       //   rho_1[ii+N*jj] = 1.0/3.0*(Rho[TV_1[i]]+Rho[TV_2[i]]+Rho[TV_3[i]]);  
       //   mu_1[ii+N*jj] = 1.0/3.0*(Mu[TV_1[i]]+Mu[TV_2[i]]+Mu[TV_3[i]]);  
       //}

       //linear interpolation
       if (dot1 >= 0 && dot2 >= 0 && dot3 >= 0) {
          //find barycentric coords
          p121 = X_11-X_31;
          p122 = X_12-X_32;
          p231 = X_21-X_31;
          p232 = X_22-X_32;

          p311 = p121*p232-p122*p231;

          chk121 = gpt1-X_31;
          chk122 = gpt2-X_32;

          chk231 = 1.0/p311*(p232*chk121-p231*chk122);
          chk232 = 1.0/p311*(p121*chk122-p122*chk121);
          chk311 = 1.0-chk231-chk232;
          
          tmpRho[k] = chk231*Rho[TV_1[i]]+chk232*Rho[TV_2[i]]+chk311*Rho[TV_3[i]];  
          tmpMu[k] = chk231*Mu[TV_1[i]]+chk232*Mu[TV_2[i]]+chk311*Mu[TV_3[i]];  
       }
    }
    }
    __kernel void calc_coefs(__global const double *X_1, __global const double *X_2, __global const unsigned *TV_1, __global const unsigned *TV_2, __global const unsigned *TV_3, __global const double *Rho, __global const double *Mu, const double dx, const unsigned N, __global double *rho_1, __global double *rho_2, __global double *mu_1, __global double *mu_2)
    { 
    unsigned i = get_global_id(0);

    //possible improvement--bring Rho, Mu, TV into local memory

    double X_11, X_12;
    double X_21, X_22;
    double X_31, X_32;

    double mn1, mn2;
    double mx1, mx2;

    double gpt1, gpt2;

    double p121, p122;
    double p231, p232;
    double p311, p312;

    double chk121, chk122;
    double chk231, chk232;
    double chk311, chk312;

    double dot1, dot2, dot3;

    int smn1, smn2; 
    int smx1, smx2; 

    uint ii,jj;

    X_11 = X_1[TV_1[i]];
    X_12 = X_2[TV_1[i]];

    X_21 = X_1[TV_2[i]];
    X_22 = X_2[TV_2[i]];

    X_31 = X_1[TV_3[i]];
    X_32 = X_2[TV_3[i]];

    //find candidate grid points
    mn1 = fmin(X_11,X_21);
    mn1 = fmin(mn1,X_31);

    mn2 = fmin(X_12,X_22);
    mn2 = fmin(mn2,X_32);

    mx1 = fmax(X_11,X_21);
    mx1 = fmax(mx1,X_31);

    mx2 = fmax(X_12,X_22);
    mx2 = fmax(mx2,X_32);

    smn1 = convert_int(ceil((mn1+0.5)/dx));
    smn2 = convert_int(ceil((mn2+0.5-dx/2.0)/dx));
    
    smx1 = convert_int(floor((mx1+0.5)/dx));
    smx2 = convert_int(floor((mx2+0.5-dx/2.0)/dx));
    
    //check whether candidate points lie inside triangle
    for (int k2 = smn2; k2 <= smx2; k2++) {
       for (int k1 = smn1; k1 <= smx1; k1++) {
          gpt1 = -0.5+k1*dx;
          gpt2 = -0.5+k2*dx+dx/2.0;

          jj = (k2 & N-1); 
          ii = (k1 & N-1); 

          //find perp vectors to half-spaces
          p121 = -(X_22-X_12);
          p122 = X_21-X_11;
          p231 = -(X_32-X_22);
          p232 = X_31-X_21;
          p311 = -(X_12-X_32);
          p312 = X_11-X_31;

          //find vectors to dot against
          chk121 = gpt1 - X_11;
          chk122 = gpt2 - X_12;
          chk231 = gpt1 - X_21;
          chk232 = gpt2 - X_22;
          chk311 = gpt1 - X_31;
          chk312 = gpt2 - X_32;

          dot1 = p121*chk121+p122*chk122;
          dot2 = p231*chk231+p232*chk232;
          dot3 = p311*chk311+p312*chk312;

          //constant interpolation
          //if (dot1 >= 0 && dot2 >= 0 && dot3 >= 0) {
          //   rho_1[ii+N*jj] = 1.0/3.0*(Rho[TV_1[i]]+Rho[TV_2[i]]+Rho[TV_3[i]]);  
          //   mu_1[ii+N*jj] = 1.0/3.0*(Mu[TV_1[i]]+Mu[TV_2[i]]+Mu[TV_3[i]]);  
          //}

          //linear interpolation
          if (dot1 >= 0 && dot2 >= 0 && dot3 >= 0) {
             //find barycentric coords
             p121 = X_11-X_31;
             p122 = X_12-X_32;
             p231 = X_21-X_31;
             p232 = X_22-X_32;

             p311 = p121*p232-p122*p231;

             chk121 = gpt1-X_31;
             chk122 = gpt2-X_32;

             chk231 = 1.0/p311*(p232*chk121-p231*chk122);
             chk232 = 1.0/p311*(p121*chk122-p122*chk121);
             chk311 = 1.0-chk231-chk232;
             
             rho_1[ii+N*jj] = chk231*Rho[TV_1[i]]+chk232*Rho[TV_2[i]]+chk311*Rho[TV_3[i]];  
             mu_1[ii+N*jj] = chk231*Mu[TV_1[i]]+chk232*Mu[TV_2[i]]+chk311*Mu[TV_3[i]];  
          }
       }
    }
    //find candidate grid points
    smn1 = convert_int(ceil((mn1+0.5-dx/2.0)/dx));
    smn2 = convert_int(ceil((mn2+0.5)/dx));
    
    smx1 = convert_int(floor((mx1+0.5-dx/2.0)/dx));
    smx2 = convert_int(floor((mx2+0.5)/dx));
    
    //check if candidate points lie inside triangle
    for (int k2 = smn2; k2 <= smx2; k2++) {
       for (int k1 = smn1; k1 <= smx1; k1++) {
          gpt1 = -0.5+k1*dx+dx/2.0;
          gpt2 = -0.5+k2*dx;

          jj = (k2 & N-1); 
          ii = (k1 & N-1); 

          //find perp vectors to half-spaces
          p121 = -(X_22-X_12);
          p122 = X_21-X_11;
          p231 = -(X_32-X_22);
          p232 = X_31-X_21;
          p311 = -(X_12-X_32);
          p312 = X_11-X_31;

          //find vectors to dot against
          chk121 = gpt1 - X_11;
          chk122 = gpt2 - X_12;
          chk231 = gpt1 - X_21;
          chk232 = gpt2 - X_22;
          chk311 = gpt1 - X_31;
          chk312 = gpt2 - X_32;

          dot1 = p121*chk121+p122*chk122;
          dot2 = p231*chk231+p232*chk232;
          dot3 = p311*chk311+p312*chk312;

          //constant interpolation
          //if (dot1 >= 0 && dot2 >= 0 && dot3 >= 0) {
          //   rho_2[ii+N*jj] = 1.0/3.0*(Rho[TV_1[i]]+Rho[TV_2[i]]+Rho[TV_3[i]]);  
          //   mu_2[ii+N*jj] = 1.0/3.0*(Mu[TV_1[i]]+Mu[TV_2[i]]+Mu[TV_3[i]]);  
          //}

          //linear interpolation
          if (dot1 >= 0 && dot2 >= 0 && dot3 >= 0) {
             //find barycentric coords
             p121 = X_11-X_31;
             p122 = X_12-X_32;
             p231 = X_21-X_31;
             p232 = X_22-X_32;

             p311 = p121*p232-p122*p231;

             chk121 = gpt1-X_31;
             chk122 = gpt2-X_32;

             chk231 = 1.0/p311*(p232*chk121-p231*chk122);
             chk232 = 1.0/p311*(p121*chk122-p122*chk121);
             chk311 = 1.0-chk231-chk232;
             
             rho_2[ii+N*jj] = chk231*Rho[TV_1[i]]+chk232*Rho[TV_2[i]]+chk311*Rho[TV_3[i]];  
             mu_2[ii+N*jj] = chk231*Mu[TV_1[i]]+chk232*Mu[TV_2[i]]+chk311*Mu[TV_3[i]];  
          }
       }
    }

    }
    __kernel void check_tri(__global const double *X_1, __global const double *X_2, __global const unsigned *TV_1, __global const unsigned *TV_2, __global const unsigned *TV_3, __global double *en)
    {
    unsigned i = get_global_id(0);

    double X_11, X_12;
    double X_21, X_22;
    double X_31, X_32;

    double a1, a2, b1, b2;
    double norm;

    X_11 = X_1[TV_1[i]];
    X_12 = X_2[TV_1[i]];

    X_21 = X_1[TV_2[i]];
    X_22 = X_2[TV_2[i]];

    X_31 = X_1[TV_3[i]];
    X_32 = X_2[TV_3[i]];

    a1 = X_21-X_11;
    a2 = X_22-X_12;

    b1 = X_31-X_11;
    b2 = X_32-X_12;

    norm = a1*b2-a2*b1;
    if (norm < 0.0) {
    en[0] = 1.0;
    }
    }
    __kernel void calc_rho(__global const double *X_1, __global const double *X_2, __global const unsigned *TV_1, __global const unsigned *TV_2, __global const unsigned *TV_3, __global const double *Rho, __global double *en)
    {
    unsigned i = get_global_id(0);

    double X_11, X_12;
    double X_21, X_22;
    double X_31, X_32;

    double a1, a2, b1, b2;
    double norm;

    X_11 = X_1[TV_1[i]];
    X_12 = X_2[TV_1[i]];

    X_21 = X_1[TV_2[i]];
    X_22 = X_2[TV_2[i]];

    X_31 = X_1[TV_3[i]];
    X_32 = X_2[TV_3[i]];

    a1 = X_21-X_11;
    a2 = X_22-X_12;

    b1 = X_31-X_11;
    b2 = X_32-X_12;

    norm = a1*b2-a2*b1;

    en[i] = 0.5*norm*1.0/3.0*(Rho[TV_1[i]]+Rho[TV_2[i]]+Rho[TV_3[i]]);
    }
    """).build()

solve_prg = cl.Program(ctx, """
    //#pragma OPENCL EXTENSION cl_amd_printf : enable
    #pragma OPENCL EXTENSION cl_khr_fp64 : enable
    //#pragma OPENCL EXTENSION cl_amd_fp64 : enable
    #define complex_ctr(x, y) (double2)(x, y)
    #define complex_mul(a, b) complex_ctr(-(a).y*(b).y+(a).x*(b).x, (a).y*(b).x+(a).x*(b).y)
    #define complex_div(a, b) 1/((b).x*(b).x+(b).y*(b).y)*complex_ctr((a).y*(b).y+(a).x*(b).x,(a).y*(b).x-(a).x*(b).y)
    __kernel void solve(__global const double2 *u_1_tf, __global const double2 *u_2_tf, __global const double2 *L_tf,
    __global const double2 *f_1_tf, __global const double2 *f_2_tf,
    const double dt, const double rho, const double mu_1, const double mu_2, const unsigned N,
    __global double2 *u_1_tf_tilde, __global double2 *u_2_tf_tilde)
    {
      unsigned i = get_global_id(0);
      unsigned j = get_global_id(1);
      unsigned base = i + N*j;
      double2 m;
      m.x = 1.0+mu_1*dt/rho*L_tf[base].x; 
      m.y = mu_1*dt/rho*L_tf[base].y; 
      double2 rhs_1 = complex_mul(m, u_1_tf[base])+dt/rho*f_1_tf[base];
      double2 rhs_2 = complex_mul(m, u_2_tf[base])+dt/rho*f_2_tf[base];
      double2 lhs;
      lhs.x = 1.0-mu_2*dt/rho*L_tf[base].x;
      lhs.y = -mu_2*dt/rho*L_tf[base].y;
      u_1_tf_tilde[base] = complex_div(rhs_1, lhs);
      u_2_tf_tilde[base] = complex_div(rhs_2, lhs);
    }
    __kernel void dal_term(__global const double *u_1, __global const double *u_2, __global const double *u_1_o, __global const double *u_2_o, __global const double *ss_1, __global const double *ss_2, __global const double *rho_1, __global const double *rho_2,
    const double dt, const unsigned N,
    __global double *f_1, __global double *f_2)
    {
      unsigned i = get_global_id(0);
      unsigned j = get_global_id(1);
      unsigned base = i + N*j;
      f_1[base] = rho_1[base]*((u_1[base]-u_1_o[base])/dt+ss_1[base]);
      f_2[base] = rho_2[base]*((u_2[base]-u_2_o[base])/dt+ss_2[base]);
    }
    __kernel void vc_laplace(__global const double *u_1,
    __global const double *u_2,
    __global const double *mu_1,
    __global const double *mu_2, const double dx, const unsigned N,
    __global double *V_1, __global double *V_2)
    {
    unsigned i = get_global_id(0);
    unsigned j = get_global_id(1);
    unsigned base = i + N*j;
    unsigned base1 = (i+1 & N-1) + N*j;
    unsigned base2 = (i-1 & N-1) + N*j;
    unsigned base3 = i + N*(j+1 & N-1);
    unsigned base4 = i + N*(j-1 & N-1);
    unsigned base5 = (i-1 & N-1) + N*(j+1 & N-1);
    unsigned base6 = (i+1 & N-1) + N*(j-1 & N-1);

    V_1[base] = 1.0/(dx*dx)*(2.0*(0.5*(mu_1[base]+mu_1[base1])*(u_1[base1]-u_1[base])
               -0.5*(mu_1[base]+mu_1[base2])*(u_1[base]-u_1[base2]))
               +0.5*(mu_1[base3]+mu_1[base])*(u_1[base3]-u_1[base])
               -0.5*(mu_1[base]+mu_1[base4])*(u_1[base]-u_1[base4]) 
               +0.5*(mu_1[base3]+mu_1[base])*(u_2[base3]-u_2[base5])
               -0.5*(mu_1[base]+mu_1[base4])*(u_2[base]-u_2[base2]));

    V_2[base] = 1.0/(dx*dx)*(0.5*(mu_2[base]+mu_2[base1])*(u_2[base1]-u_2[base])
               -0.5*(mu_2[base]+mu_2[base2])*(u_2[base]-u_2[base2])
               +0.5*(mu_2[base]+mu_2[base1])*(u_1[base1]-u_1[base6])
               -0.5*(mu_2[base]+mu_2[base2])*(u_1[base]-u_1[base4])
               +2*(0.5*(mu_2[base]+mu_2[base3])*(u_2[base3]-u_2[base])
               -0.5*(mu_2[base]+mu_2[base4])*(u_2[base]-u_2[base4])));

    }   
    """).build()

prg = cl.Program(ctx, """
    #pragma OPENCL EXTENSION cl_khr_fp64 : enable
    //#pragma OPENCL EXTENSION cl_amd_fp64 : enable
    __kernel void p_filt(__global double2 *u)
    {
      u[0].x = 0;
      u[0].y = 0;
    }
    __kernel void L_adj_1(__global double2 *L_tf)
    {
      L_tf[0].x = 0.0;
      L_tf[0].y = 0.0;
    }
    __kernel void L_adj_2(__global double2 *L_tf)
    {
      L_tf[0].x = 1.0;
      L_tf[0].y = 0.0;
    }
    __kernel void nyq_filt(__global double2 *u, const unsigned N)
    {
      unsigned i = get_global_id(0);
      unsigned base = N/2 + N*i; 
      unsigned base1 = i + N*N/2; 
   
      u[base].x = 0;
      u[base].y = 0;
      u[base1].x = 0;
      u[base1].y = 0;
    }
    __kernel void vort(__global const double *u_1,
    __global const double *u_2,
    const double dx, const unsigned N, __global double *f)
    {
      unsigned i = get_global_id(0);
      unsigned j = get_global_id(1);
      unsigned base = i + N*j;
      unsigned base1 = (i-1 & N-1)+N*j;
      unsigned base2 = i+N*(j-1 & N-1);
      f[base] = 1/dx*(u_2[base]-u_2[base1]-(u_1[base]-u_1[base2]));
    }
    __kernel void div(__global const double *u_1,
    __global const double *u_2,
    const double dx, const unsigned N, __global double *f)
    {
      unsigned i = get_global_id(0);
      unsigned j = get_global_id(1);
      //f[i,j] = 1/dx*(u_1[i_2[i,j],ip_1[i,j]]-u_1[i,j]+u_2[ip_2[i,j],i_1[i,j]]-u_2[i,j];
      unsigned base = i + N*j;
      unsigned base1 = (i+1 & N-1)+N*j;
      unsigned base2 = i+N*(j+1 & N-1);
      f[base] = 1/dx*(u_1[base1]-u_1[base]+u_2[base2]-u_2[base]);
    }
    __kernel void grad(__global const double *phi,
    const double dx, const unsigned N,
    __global double *g_1, __global double *g_2)
    {
      unsigned i = get_global_id(0);
      unsigned j = get_global_id(1);
      //grad_1 = 1/dx*(phi[i_2,i_1]-phi[i_2,im_1])
      //grad_2 = 1/dx*(phi[i_2,i_1]-phi[im_2,i_1])
      unsigned base = i+N*j;
      unsigned base1 = (i-1 & N-1)+N*j;
      unsigned base2 = i+N*(j-1 & N-1);
      g_1[base] = 1/dx*(phi[base]-phi[base1]);
      g_2[base] = 1/dx*(phi[base]-phi[base2]);
    }
    __kernel void ss_conv(__global const double *u_1,
    __global const double *u_2, const double dx, const unsigned N,
    __global double *S_1, __global double *S_2)
    {
      unsigned i = get_global_id(0);
      unsigned j = get_global_id(1);
      //Compute S = 0.5*(div(UU) + (U*grad)U).
      unsigned base = i+N*j;
      unsigned base1 = (i-1 & N-1)+N*j;
      unsigned base2 = (i+1 & N-1)+N*j;
      unsigned base3 = i+N*(j-1 & N-1);
      unsigned base4 = i+N*(j+1 & N-1);
      unsigned base5 = (i-1 & N-1)+N*(j+1 & N-1);
      unsigned base6 = (i+1 & N-1)+N*(j-1 & N-1);
      double D_x, D_y, A_x, A_y, D, A;
      D_x = (0.25/dx)*(
          (u_1[base2]+u_1[base])*(u_1[base2]+u_1[base]) -
          (u_1[base]+u_1[base1])*(u_1[base]+u_1[base1]));
      D_y = (0.25/dx)*(
          (u_2[base4]+u_2[base5])*(u_1[base4]+u_1[base]) -
          (u_2[base]+u_2[base1])*(u_1[base]+u_1[base3]));
      D = D_x + D_y;
      A_x =  (0.25/dx)*(
          (u_1[base2]+u_1[base])*(u_1[base2]-u_1[base]) +
          (u_1[base]+u_1[base1])*(u_1[base]-u_1[base1]));
      A_y = (0.25/dx)*(
            (u_2[base4]+u_2[base5])*(u_1[base4]-u_1[base]) +
            (u_2[base]+u_2[base1])*(u_1[base]-u_1[base3]));
      A = A_x + A_y;
      S_1[base] = 0.5*(D+A);
      D_x = (0.25/dx)*(
          (u_1[base2]+u_1[base6])*(u_2[base2]+u_2[base]) -
          (u_1[base]+u_1[base3])*(u_2[base]+u_2[base1]));
      D_y = (0.25/dx)*(
          (u_2[base4]+u_2[base])*(u_2[base4]+u_2[base]) -
          (u_2[base]+u_2[base3])*(u_2[base]+u_2[base3]));
      D = D_x + D_y;
      A_x =  (0.25/dx)*(
          (u_1[base2]+u_1[base6])*(u_2[base2]-u_2[base]) +
          (u_1[base]+u_1[base3])*(u_2[base]-u_2[base1]));
      A_y = (0.25/dx)*(
            (u_2[base4]+u_2[base])*(u_2[base4]-u_2[base]) +
            (u_2[base]+u_2[base3])*(u_2[base]-u_2[base3]));
      A = A_x + A_y;
      S_2[base] = 0.5*(D+A);
    }
    __kernel void interpctolr(__global const double *phi, __global const unsigned *i_2,  __global const unsigned *im_1, const unsigned N, __global double *f)
    {
      unsigned i = get_global_id(0);
      unsigned j = get_global_id(1);
      //phi_itp = 0.5*(phi+phi[i_2,im_1])

      unsigned base = i+N*j;
      unsigned base1 = im_1[base]+N*i_2[base];
      f[base] = 0.5*(phi[base]+phi[base1]);
    }
    __kernel void interpctotb(__global const double *phi, __global const unsigned *i_1, __global const unsigned *im_2, const unsigned N, __global double *f)
    {
      unsigned i = get_global_id(0);
      unsigned j = get_global_id(1);
      //phi_itp = 0.5*(phi+phi[im_2,i_1])

      unsigned base = i+N*j;
      unsigned base1 = i_1[base]+N*im_2[base];
      f[base] = 0.5*(phi[base]+phi[base1]);
    }
    __kernel void interplrtoc(__global const double *phi, __global const unsigned *i_2, __global const unsigned *ip_1, const unsigned N, __global double *f)
    {
      unsigned i = get_global_id(0);
      unsigned j = get_global_id(1);
      //phi_itp = 0.5*(phi+phi[i_2,ip_1])

      unsigned base = i+N*j;
      unsigned base1 = ip_1[base]+N*i_2[base];
      f[base] = 0.5*(phi[base]+phi[base1]);
    }
    __kernel void interptbtoc(__global const double *phi, __global const unsigned *i_1,__global const unsigned *ip_2, const unsigned N, __global double *f)
    {
      unsigned i = get_global_id(0);
      unsigned j = get_global_id(1);
      //phi_itp = 0.5*(phi+phi[ip_2,i_1])

      unsigned base = i+N*j;
      unsigned base1 = i_1[base]+N*ip_2[base];
      f[base] = 0.5*(phi[base]+phi[base1]);
    }
    """).build()

calc_lagF_teth_kernel = ib_prg.calc_lagF_teth
calc_lagF_teth_kernel.set_scalar_arg_dtypes([None]*4 + [np.uint32] + [np.float64]*2 + [None]*2)

calc_lagF_kernel = ib_prg.calc_lagF
calc_lagF_kernel.set_scalar_arg_dtypes([None]*5 + [np.float64]*3 + [None]*2)

div_kernel = prg.div
div_kernel.set_scalar_arg_dtypes([None]*2 + [np.float64, np.uint32] + [None])

vort_kernel = prg.vort
vort_kernel.set_scalar_arg_dtypes([None]*2 + [np.float64, np.uint32] + [None])

L_adj_1_kernel = prg.L_adj_1
L_adj_1_kernel.set_scalar_arg_dtypes([None])

L_adj_2_kernel = prg.L_adj_2
L_adj_2_kernel.set_scalar_arg_dtypes([None])

p_filt_kernel = prg.p_filt
p_filt_kernel.set_scalar_arg_dtypes([None])

nyq_filt_kernel = prg.nyq_filt
nyq_filt_kernel.set_scalar_arg_dtypes([None] + [np.uint32])

solve_kernel = solve_prg.solve
solve_kernel.set_scalar_arg_dtypes([None]*5 + [np.float64]*4 + [np.uint32] + [None]*2)

vc_laplace_kernel = solve_prg.vc_laplace
vc_laplace_kernel.set_scalar_arg_dtypes([None]*4 + [np.float64] + [np.uint32] + [None]*2)

dal_term_kernel = solve_prg.dal_term
dal_term_kernel.set_scalar_arg_dtypes([None]*8 + [np.float64] + [np.uint32] + [None]*2)

grad_kernel = prg.grad
grad_kernel.set_scalar_arg_dtypes([None]*1 + [np.float64, np.uint32] + [None]*2)

ss_conv_kernel = prg.ss_conv
ss_conv_kernel.set_scalar_arg_dtypes([None]*2 + [np.float64, np.uint32] + [None]*2)

interp_kernel = ib_prg.interp
interp_kernel.set_scalar_arg_dtypes([None]*4 + [np.uint32, np.float64] + [None]*2)
spread_kernel = ib_prg.spread
spread_kernel.set_scalar_arg_dtypes([None]*8 + [np.uint32] + [np.float64]*2 + [np.uint32]*2 + [None]*2)
updaten_kernel = ib_prg.updaten
updaten_kernel.set_scalar_arg_dtypes([None]*6 + [np.uint32] + [np.float64] + [np.uint32]*2)
calc_coefs_retriang_kernel = ib_prg.calc_coefs_retriang
calc_coefs_retriang_kernel.set_scalar_arg_dtypes([None]*7 + [np.uint32]*2 + [None]*2)
calc_coefs_kernel = ib_prg.calc_coefs
calc_coefs_kernel.set_scalar_arg_dtypes([None]*7 + [np.float64] + [np.uint32] + [None]*4)
calc_rho_kernel = ib_prg.calc_rho
calc_rho_kernel.set_scalar_arg_dtypes([None]*7)

check_tri_kernel = ib_prg.check_tri
check_tri_kernel.set_scalar_arg_dtypes([None]*6)

interpctolr_kernel = prg.interpctolr
interpctolr_kernel.set_scalar_arg_dtypes([None]*3 + [np.uint32] + [None])
interpctotb_kernel = prg.interpctotb
interpctotb_kernel.set_scalar_arg_dtypes([None]*3 + [np.uint32] + [None])
interplrtoc_kernel = prg.interplrtoc
interplrtoc_kernel.set_scalar_arg_dtypes([None]*3 + [np.uint32] + [None])
interptbtoc_kernel = prg.interptbtoc
interptbtoc_kernel.set_scalar_arg_dtypes([None]*3 + [np.uint32] + [None])

for ii in range(start_size,end_size+1): 
    #set up grid
    N = 2**ii
    dx = 1.0/N
    rho_0 = 0.0

    x_1_stag = np.arange(-0.5, 0.5, dx)
    x_2_stag = np.arange(-0.5, 0.5, dx)
    xx_1_stag, xx_2_stag = np.meshgrid(x_1_stag, x_2_stag)

    x_1_c = np.arange(-0.5+dx/2, 0.5+dx/2, dx)
    x_2_c = np.arange(-0.5+dx/2, 0.5+dx/2, dx)
    xx_1_c, xx_2_c = np.meshgrid(x_1_c, x_2_c)
    
    #region = np.zeros((N,N))
    region = np.zeros((N*N,1))
    y_cen = 0.1
    y_cen = 0.15
    #y_cen = 0.2
    y_cen = 0.0
    #y_width = 0.5
    y_width = 0.05
    #y_width = 0.15
    #y_width = 0.25
    #y_width = 0.3
    #y_width = 0.2
    #y_width = 0.25
    #y_width = 0.4
    #y_width = 0.1
    #y_width = 0.3
    #y_width = 0.01
    #y_width = 0.1
    #y_width = 0.2
    x_cen = 0.0
    #x_cen = -0.2
    #x_cen = -0.1
    #x_cen = -0.15
    x_width = 0.4
    x_width = 0.5
    x_width = 0.05
    #x_width = 0.3
    #x_width = 0.2
    #x_width = 0.2
    #x_width = 0.15
    #x_width = 0.25
    #x_width = 0.3
    #x_width = 0.5
    #x_width = 0.4
    #x_width = 0.1
    #x_width = 0.2
    #inds = np.nonzero(np.logical_or(xx_2_c > y_cen+y_width,xx_2_c < y_cen-y_width))
    #for i in range(0,len(inds[0])):
    #    for j in range(0,len(inds[1])):
    #        region[inds[0][i],inds[1][j]] = 0.5*(1.0+np.cos(pi*(xx_2_c[inds[0][i],inds[1][j]]-y_cen)/y_width))
    xx_2_c_flat = xx_2_c.reshape(N*N,1)
    xx_1_stag_flat = xx_1_stag.reshape(N*N,1)
    region[(np.logical_and(np.logical_and(xx_1_stag_flat < x_cen+x_width,xx_1_stag_flat > x_cen-x_width),np.logical_and(xx_2_c_flat < y_cen+y_width,xx_2_c_flat > y_cen-y_width)))] = 0.5*(1.0+np.cos(pi*(xx_2_c_flat[(np.logical_and(np.logical_and(xx_1_stag_flat < x_cen+x_width,xx_1_stag_flat > x_cen-x_width),np.logical_and(xx_2_c_flat < y_cen+y_width,xx_2_c_flat > y_cen-y_width)))]-y_cen)/y_width))*0.5*(1.0+np.cos(pi*(xx_1_stag_flat[(np.logical_and(np.logical_and(xx_1_stag_flat < x_cen+x_width,xx_1_stag_flat > x_cen-x_width),np.logical_and(xx_2_c_flat < y_cen+y_width,xx_2_c_flat > y_cen-y_width)))]-x_cen)/x_width))
    #region[(np.logical_and(np.logical_and(xx_1_stag_flat < x_cen+x_width,xx_1_stag_flat > x_cen-x_width),np.logical_and(xx_2_c_flat < y_cen+y_width,xx_2_c_flat > y_cen-y_width)))] = 1.0 
    #region = 1.0/(x_width*y_width*2.0*pi)*np.exp(-(xx_2_c_flat-y_cen)**2/(2*y_width**2))*np.exp(-(xx_1_stag_flat-x_cen)**2/(2*x_width**2))
    #region[(np.logical_and(np.logical_and(xx_1_stag_flat < x_cen+x_width,xx_1_stag_flat > x_cen-x_width),np.logical_and(xx_2_c_flat < y_cen+y_width,xx_2_c_flat > y_cen-y_width)))] = 1.0 
    #region[(np.logical_and(xx_2_c_flat < y_cen+y_width,xx_2_c_flat > y_cen-y_width))] = 0.5*(1.0+np.cos(pi*(xx_2_c_flat[(np.logical_and(xx_2_c_flat < y_cen+y_width,xx_2_c_flat > y_cen-y_width))]-y_cen)/y_width))
    #region[(np.logical_and(xx_2_c_flat < y_cen+y_width,xx_2_c_flat > y_cen-y_width))] = 1.0 
    sqregion = np.zeros((N*N,1))
    #sqregion[(np.logical_and(xx_2_c_flat < y_cen+y_width,xx_2_c_flat > y_cen-y_width))] = 1.0 
    #sqregion[(np.logical_and(np.logical_and(xx_1_stag_flat < x_cen+x_width,xx_1_stag_flat > x_cen-x_width),np.logical_and(xx_2_c_flat < y_cen+y_width,xx_2_c_flat > y_cen-y_width)))] = 1.0 
    #sqregion[(np.logical_and(xx_2_c_flat < 0.45,xx_2_c_flat > -0.45))] = 1.0 
    sqregion[(np.logical_and(np.logical_and(xx_1_stag_flat < 0.4,xx_1_stag_flat > -0.4),np.logical_and(xx_2_c_flat < 0.4,xx_2_c_flat > -0.4)))] = 1.0 
    #sqregion[(np.logical_and(xx_2_c_flat < 0.5-dx,xx_2_c_flat > -0.5+dx))] = 1.0 
    sqregion = np.ones((N*N,1))
    sqregion = sqregion.reshape(N,N)
    #region = np.ones((N*N,1))
    #region = np.ones((N*N,1))
    region = region.reshape(N,N)
    #inds = np.nonzero(1-np.logical_or(xx_1_c > -0.0,np.logical_or(np.logical_or(xx_1_c < -0.4,xx_2_c > 0.4),xx_2_c < 0.0)))
    #inds = np.nonzero(1-np.logical_or(xx_1_c > -0.0,np.logical_or(np.logical_or(xx_1_c < -0.4,xx_2_c > 0.4),xx_2_c < 0.0)))
    #region[np.nonzero(np.logical_or(xx_1_c > -0.1,np.logical_or(np.logical_or(xx_1_c < -0.3,xx_2_c > 0.4),xx_2_c < 0.2)))] = 0
    #region[np.logical_or(xx_1_c > -0.0,np.logical_or(np.logical_or(xx_1_c < -0.4,xx_2_c > 0.4),xx_2_c < 0.0)] = 0
    #region[np.ix_(inds[0],inds[1])] = 1
    #region[inds] = 0

    if PR_INT:
        mesh3 = np.vstack((xx_1_c[np.newaxis],xx_2_c[np.newaxis],np.zeros((N,N))[np.newaxis]))

    mu_0 = 0.05
    mu_0 = 0.01
    mu_0 = 0.001
    #mu_0 = 0.0001
    mu_0 = 0.00001
    #mu_0 = 0.000001
    #mu_0 = 0.0000001
    mu_0 = 0.000000001
    #mu_0 = 0.00000000001
    #mu_0 = 0.0000000000000001
    #mu_0 = 0.0
    #mu_0 = 0.005
    dt = 0.1*dx
    #dt = 0.025*dx
    #dt = 0.01*dx
    #dt = 0.05*dx
    clock_max = np.int(np.round(0.4/dt,0))
    #clock_max = np.int(np.round(10.0/dt,0))
    clock_max = np.int(np.round(100.0/dt,0))
    #clock_max = np.int(np.round(60.0/dt,0))
    #clock_max = np.int(np.round(20.0/dt,0))
    #clock_max = np.int(np.round(16.0/dt,0))
    #clock_max = np.int(np.round(4.0/dt,0))
    #clock_max = 2
    print 'grid size', N
    print 'dt is', dt
    print 'total timesteps', clock_max

    #set up Lagrangian parameters
    myratio = 3./4
    #N_r = N-1
    #N_s = N-1
    N_r = N/2-1
    N_s = N/2-1
    #N_r = N/2
    #N_s = N/2
    #N_r = N-1
    #N_s = N-1
    #N_r = 2*N-1
    #N_s = 2*N-1
    #N_r = N/4-1
    #N_s = N/4-1
    #N_r = N/4-1
    #N_s = N/4-1
    #N_r = N/2-2
    #N_s = N/2-2
    #N_r = N/2
    #N_s = N/2
    #N_r = N/4
    #N_s = N/4
    #dr = 1./N_r
    #ds = 1./N_s
    #debugging
    #dr = (1.0-0.00001)/N_r
    #ds = (1.0-0.00001)/N_s
    dr = 1.0/N_r
    ds = 1.0/N_s
    dtheta = dr*ds
    N_r = N_r+1
    N_s = N_s+1
    #N_r = N_r+2
    #N_s = N_s+2
    print 'N_r is', N_r
    b_q = N_r*N_s
    b_add = 0
    N_edge = 0

    r_lag = np.arange(-0.5+dr/2, 0.5+dr/2, dr)
    s_lag = np.arange(-0.5+ds/2, 0.5+ds/2, ds)
    #r_lag = np.arange(dr/2, 1.0+dr/2, dr)
    #s_lag = np.arange(ds/2, 1.0+ds/2, ds)
    r_lag = np.arange(dr/2, 1.0+3*dr/2, dr)
    s_lag = np.arange(ds/2, 1.0+3*ds/2, ds)
    r_lag = np.arange(-0.000005, (1.0-0.000005)+dr, dr)
    s_lag = np.arange(-0.000005, (1.0-0.000005)+ds, ds)
    #r_lag = np.arange(-dr/2, 1.0+dr/2, dr)
    #s_lag = np.arange(-ds/2, 1.0+ds/2, ds)
    rr_lag, ss_lag = np.meshgrid(r_lag, s_lag)
    #ss_lag, rr_lag = np.meshgrid(s_lag, r_lag)
    rr_lag = rr_lag.reshape(b_q,1).flatten(1)
    ss_lag = ss_lag.reshape(b_q,1).flatten(1)
    #rr_lag = (rr_lag.T).reshape(b_q,1).flatten(1)
    #ss_lag = (ss_lag.T).reshape(b_q,1).flatten(1)

    X_1 = (rr_lag-0.5).copy()
    X_2 = (ss_lag-0.5).copy()
    #X_1 = (myratio*rr_lag-0.5).copy()
    #X_2 = (myratio*ss_lag-0.5).copy()
    X_1 = X_1.reshape(b_q,1).copy()
    X_2 = X_2.reshape(b_q,1).copy()

    lagregion = np.zeros((b_q,1))
    #lagregion[X_1 > -0.1 or X_1 < -0.3 or X_2 > 0.4 or X_2 < 0.2] = 0
    #lagregion[np.nonzero(np.logical_or(X_1 > -0.1,np.logical_or(np.logical_or(X_1 < -0.3,X_2 > 0.4),X_2 < 0.2))] = 0
    #inds = np.nonzero(1-np.logical_or(X_1 > -0.0,np.logical_or(np.logical_or(X_1 < -0.4,X_2 > 0.4),X_2 < 0.0)))
    #lagregion[inds[0],inds[1]] = 0
    #lagregion[inds] = 0
    #lagregion[np.nonzero(np.logical_or(X_1 > -0.1,np.logical_or(np.logical_or(X_1 < -0.3,X_2 > 0.4),X_2 < 0.2)))] = 0
    #lagregion[inds[0]] = 1
    #inds = np.nonzero(np.logical_or(X_2 > y_cen+y_width,X_2 < y_cen-y_width))
    #for i in range(0,len(inds[0])):
    #    region[inds[0][i]] = 1.0+np.cos(pi*(X_2[inds[0][i]]-y_cen)/y_width)
    #lagregion[(np.logical_and(X_2 < y_cen+y_width,X_2 > y_cen-y_width))] = 0.5*(1.0+np.cos(pi*(X_2[(np.logical_and(X_2 < y_cen+y_width,X_2 > y_cen-y_width))]-y_cen)/y_width))
    lagregion[(np.logical_and(np.logical_and(X_1 < x_cen+x_width,X_1 > x_cen-x_width),np.logical_and(X_2 < y_cen+y_width,X_2 > y_cen-y_width)))] = 0.5*(1.0+np.cos(pi*(X_2[(np.logical_and(np.logical_and(X_1 < x_cen+x_width,X_1 > x_cen-x_width),np.logical_and(X_2 < y_cen+y_width,X_2 > y_cen-y_width)))]-y_cen)/y_width))*0.5*(1.0+np.cos(pi*(X_1[(np.logical_and(np.logical_and(X_1 < x_cen+x_width,X_1 > x_cen-x_width),np.logical_and(X_2 < y_cen+y_width,X_2 > y_cen-y_width)))]-x_cen)/x_width))
    #lagregion[(np.logical_and(np.logical_and(X_1 < x_cen+x_width,X_1 > x_cen-x_width),np.logical_and(X_2 < y_cen+y_width,X_2 > y_cen-y_width)))] = 1.0 
    #lagregion = 1.0/(x_width*y_width*2.0*pi)*np.exp(-(X_2-y_cen)**2/(2*y_width**2))*np.exp(-(X_1-x_cen)**2/(2*x_width**2))
    #lagregion[(np.logical_and(np.logical_and(X_1 < x_cen+x_width,X_1 > x_cen-x_width),np.logical_and(X_2 < y_cen+y_width,X_2 > y_cen-y_width)))] = 1.0 
    #lagregion[(np.logical_and(X_2 < y_cen+y_width,X_2 > y_cen-y_width))] = 1.0 
    #lagregion = np.ones((b_q,1))
    #lagregion = np.ones((b_q,1))
    lagregion = lagregion.flatten(1)
    #lagregion = (lagregion.T).flatten(1)

    N_tri = 2*(N_r-1)*(N_s-1)
    #N_tri = 2*(int(np.sqrt(b_q))-1)*int(np.sqrt(b_q))
    TV = np.zeros((N_tri,3))

    for i in range(0,N_r-1):
        for j in range(0,N_s-1):
            #TV[2*(j+(np.sqrt(b_q)-1)*i),0] = j+np.sqrt(b_q)*i
            #TV[2*(j+(np.sqrt(b_q)-1)*i),2] = (np.mod(np.mod(j+np.sqrt(b_q)*i+1,b_q),np.sqrt(b_q))+np.sqrt(b_q)*np.floor(np.mod(j+np.sqrt(b_q)*i,b_q)/np.sqrt(b_q)))
            #TV[2*(j+(np.sqrt(b_q)-1)*i),1] = np.mod(j+np.sqrt(b_q)*i+np.sqrt(b_q),b_q)
            #TV[2*(j+(np.sqrt(b_q)-1)*i)+1,0] = np.mod(j+np.sqrt(b_q)*i+np.sqrt(b_q),b_q)
            #TV[2*(j+(np.sqrt(b_q)-1)*i)+1,2] = (np.mod(np.mod(j+np.sqrt(b_q)*i+1,b_q),np.sqrt(b_q))+np.sqrt(b_q)*np.floor(np.mod(j+np.sqrt(b_q)*i,b_q)/np.sqrt(b_q)))
            #TV[2*(j+(np.sqrt(b_q)-1)*i)+1,1] = np.mod(np.mod(np.mod(j+np.sqrt(b_q)*i+1,b_q),np.sqrt(b_q))+np.sqrt(b_q)*np.floor(np.mod(j+np.sqrt(b_q)*i,b_q)/np.sqrt(b_q))+np.sqrt(b_q),b_q)
            #TV[2*(j+(N_s-1)*i),0] = j+N_s*i
            #TV[2*(j+(N_s-1)*i),2] = j+N_s*i+1
            #TV[2*(j+(N_s-1)*i),1] = j+N_s*i+N_s
            #TV[2*(j+(N_s-1)*i)+1,0] = j+N_s*i+N_s
            #TV[2*(j+(N_s-1)*i)+1,2] = j+N_s*i+1
            #TV[2*(j+(N_s-1)*i)+1,1] = j+N_s*i+1+N_s
            TV[2*(i+(N_r-1)*j),0] = i+N_r*j
            TV[2*(i+(N_r-1)*j),1] = i+N_r*j+1
            TV[2*(i+(N_r-1)*j),2] = i+N_r*j+N_r
            TV[2*(i+(N_r-1)*j)+1,0] = i+N_r*j+N_r
            TV[2*(i+(N_r-1)*j)+1,1] = i+N_r*j+1
            TV[2*(i+(N_r-1)*j)+1,2] = i+N_r*j+1+N_r
            #TV[2*(j+(N_s-1)*i),0] = j+N_s*i
            #TV[2*(j+(N_s-1)*i),2] = j+N_s*i+1
            #TV[2*(j+(N_s-1)*i),1] = j+N_s*i+N_s
            #TV[2*(j+(N_s-1)*i)+1,0] = j+N_s*i+N_s
            #TV[2*(j+(N_s-1)*i)+1,2] = j+N_s*i+1
            #TV[2*(j+(N_s-1)*i)+1,1] = j+N_s*i+1+N_s

    edges = []
    for i in range(0, N_s-1):
      edges.append((i, i+1))
    for i in range(0, N_s-1):
      edges.append((N_s*(N_r-1)+i, N_s*(N_r-1)+i+1))
    for i in range(0, N_r-1):
      edges.append((N_s*i, N_s*i+N_s))
    for i in range(0, N_r-1):
      edges.append((N_s*(i+1)-1, N_s*(i+1)+N_s-1))
    edges_array = np.asarray(edges)

    mu_l = 10.0
    mu_l = 1.0
    mu_l = 0.1
    mu_l = 0.0
    rho_l = 10.0
    rho_l = 0.1
    rho_l = 1.0
    #rho_l = 2.0
    #rho_l = 0.0
    #Mu = mu_l*0.5*mu_0*(1.+np.sin(2*pi*rr_lag-pi/2))
    #Mu = mu_0*np.ones((b_q,1))
    #Rho = rho_l*0.5*rho_0*(1.+np.sin(2*pi*rr_lag-pi/2))
    #Rho = rho_0*np.ones((b_q,1))
    grav = 1000.0 
    grav = 0.000001 
    grav = 0.0001
    grav = 0.001
    grav = 1.0
    grav = 10.0
    grav = 2.0
    #grav = 200.0
    #grav = 40.0
    #grav = 0.1
    grav = 100.0
    grav = 400.0
    #grav = 4000.0
    #grav = 40000.0
    #grav = 40.0
    #grav = 4.0
    #grav = 1000.0
    #grav = 1000.0
    #grav = 400.0
    #grav = 0.005
    #grav = 10.0 
    body_f = 0.01
    #body_f = 0.025
    #body_f = 0.2
    #body_f = 1.0
    body_f = 1.0
    #body_f = 0.001
    #body_f = 0.0001
    #body_f = 0.0

    C_v = 0.01
    C_v = 0.001
    #C_v = 0.0001
    #C_v = 0.000001
    #C_v = 0.0001
    #C_v = 0.01
    #C_v = 0.00001
    #C_v = 0.01
    C_v = 0.0
    C_h = 12*pi
    C_h = 48*pi
    C_h = grav
    C_h = 24
    #C_h = 48
    #C_h = 480*pi
    #C_h = 2*pi
    #C_h = 4.8*pi
    #C_h = 0.48*pi
    #wave_k = 4.0
    wave_l = 3.0
    wave_k = 3.0
    wave_l = 6.0
    wave_k = 6.0
    wave_l = 1.0
    wave_l = 40.0
    wave_l = 6.0
    #wave_l = 0.0
    #wave_l = 32.0
    wave_k = 20.0
    wave_k = 40.0
    wave_k = 1.0
    wave_k = 6.0
    #wave_k = 32.0
    #wave_l = 32.0
    #wave_k = 32.0
    #wave_k = 32.0
    #wave_l = 3.0
    #wave_l = 16.0
    #wave_l = 32.0
    #wave_k = 1.0
    #wave_l = 0.0
    #buoy_freq = np.sqrt(grav/(2*pi*C_h))
    buoy_freq = np.sqrt(grav/C_h)
    print 'buoy freq is', buoy_freq
    #buoy_freq = np.sqrt(grav/C_h)
    #buoy_freq = np.sqrt(grav/C_h)
    wave_omega = buoy_freq*wave_k/np.sqrt(wave_k*wave_k+wave_l*wave_l+1.0/(4*C_h*C_h)) 
    theta = pi/3.0
    theta = pi/4.0
    wave_omega = buoy_freq*np.cos(theta)
    #print 'period is', 1.0/wave_omega 
    print 'period is', 2.0*pi/wave_omega 
    u_2_soln = C_v*np.exp(2*pi*xx_2_stag/(2*C_h))*(np.exp(1j*2*pi*(wave_k*xx_1_c+wave_l*xx_2_stag))).real*region
    #u_2_soln = C_v*np.exp(xx_2_stag/(2*C_h))*(np.exp(1j*(wave_k*xx_1_c+wave_l*xx_2_stag))).real*region
    #u_1_soln = ((2.0*pi*1j)/wave_k*(1.0/(2*C_h)+1j*wave_l)*C_v*np.exp(2*pi*xx_2_c/(2*C_h))*np.exp(1j*(2*pi*wave_k*xx_1_stag+2*pi*wave_l*xx_2_c))).real*region
    u_1_soln = (1j/wave_k*(1j*wave_l+1.0/(2*C_h))*C_v*np.exp(2*pi*xx_2_c/(2*C_h))*np.exp(1j*2*pi*(wave_k*xx_1_stag+wave_l*xx_2_c))).real*region
    #u_1_soln = (1j/wave_k*(1j*wave_l+1.0/(2*C_h))*C_v*np.exp(xx_2_c/(2*C_h))*np.exp(1j*(wave_k*xx_1_stag+wave_l*xx_2_c))).real*region
    #u_2_old_soln = C_v*np.exp(2*pi*xx_2_stag/(2*C_h))*(np.exp(1j*2*pi*(wave_k*xx_1_c+wave_l*xx_2_stag+wave_omega*dt))).real*region
    #u_1_old_soln = (1j/wave_k*(1j*wave_l+1.0/(2*C_h))*C_v*np.exp(2*pi*xx_2_c/(2*C_h))*np.exp(1j*2*pi*(wave_k*xx_1_stag+wave_l*xx_2_c+wave_omega*dt))).real*region
    u_2_old_soln = C_v*np.exp(2*pi*xx_2_stag/(2*C_h))*(np.exp(1j*2*pi*(wave_k*xx_1_c+wave_l*xx_2_stag+wave_omega*dt))).real*region
    u_1_old_soln = (1j/wave_k*(1j*wave_l+1.0/(2*C_h))*C_v*np.exp(2*pi*xx_2_c/(2*C_h))*np.exp(1j*2*pi*(wave_k*xx_1_stag+wave_l*xx_2_c+wave_omega*dt))).real*region
    #u_2_old_soln = C_v*np.exp(xx_2_stag/(2*C_h))*(np.exp(1j*(wave_k*xx_1_c+wave_l*xx_2_stag+wave_omega*dt))).real*region
    #u_1_old_soln = (1j/wave_k*(1j*wave_l+1.0/(2*C_h))*C_v*np.exp(xx_2_c/(2*C_h))*np.exp(1j*(wave_k*xx_1_stag+wave_l*xx_2_c+wave_omega*dt))).real*region

    #Mu = mu_l*0.5*(1.+np.sin(2*pi*ss_lag-pi))*0.5*mu_0*(1.+np.sin(2*pi*rr_lag))
    #Rho = rho_l*0.5*(1.+np.sin(2*pi*ss_lag-pi))*0.5*rho_0*(1.+np.sin(2*pi*rr_lag))
    #Mu = mu_l*0.5*(1.+np.sin(2*pi*ss_lag-pi)*np.sin(2*pi*rr_lag))*mu_0
    #Rho = rho_l*0.5*(1.+np.sin(2*pi*ss_lag-pi)*np.sin(2*pi*rr_lag))*rho_0
    Mu = mu_l*0.25*((1.0-ss_lag)*(1.0+0.5*np.sin(pi*(rr_lag-0.1*ss_lag))))*mu_0
    Rho = rho_l*0.25*((1.0-ss_lag)*(1.0+0.5*np.sin(pi*(rr_lag-0.1*ss_lag))))*rho_0
    #Mu = mu_l*0.5*((1.0-ss_lag)*(1.0+(rr_lag-0.1*ss_lag)))*mu_0
    #Rho = rho_l*0.5*((1.0-ss_lag)*(1.0+(rr_lag-0.1*ss_lag)))*rho_0
    Mu = mu_l*1./6.*((1.0-ss_lag)*(1.0+5.0*(rr_lag-0.1*ss_lag)))*mu_0
    Rho = rho_l*1./6.*((1.0-ss_lag)*(1.0+5.0*(rr_lag-0.1*ss_lag)))*rho_0
    #Mu = mu_l*0.25*(1.+np.sin(2*pi*ss_lag-pi/2.))*(1.+np.sin(2*pi*rr_lag-pi/2.))*mu_0
    #Rho = rho_l*0.25*(1.+np.sin(2*pi*ss_lag-pi/2.))*(1.+np.sin(2*pi*rr_lag-pi/2.))*rho_0
    Rho = rho_l*rho_0*np.exp(-ss_lag)+0.1*rho_l*rho_0*np.exp(0.5*(-ss_lag))*np.sin(4*pi*(-rr_lag)-pi/2)*np.sin(2*pi*(-ss_lag)-pi/2)
    Mu = mu_l*mu_0*np.exp(-ss_lag)+0.1*mu_l*mu_0*np.exp(0.5*(-ss_lag))*np.sin(4*pi*(-rr_lag)-pi/2)*np.sin(2*pi*(-ss_lag)-pi/2)
    Rho = rho_l*rho_0*np.exp(-ss_lag)+1.0*rho_l*rho_0*np.exp(-50.0*((rr_lag-0.25)**2+(ss_lag-0.75)**2))
    Mu = mu_l*mu_0*np.exp(-ss_lag)+1.0*mu_l*mu_0*np.exp(-50.0*((rr_lag-0.25)**2+(ss_lag-0.75)**2))
    Mu = mu_l*mu_0*np.exp(-ss_lag)+0.1*mu_l*mu_0*(1./(2.*pi*1./50))**0.25*np.cos(2.*pi*(rr_lag-0.25)+2.*pi*(ss_lag-0.75))*np.exp(0.25/(1./50)*-((rr_lag-0.25)**2+(ss_lag-0.75)**2))
    Rho = rho_l*rho_0*np.exp(-ss_lag)+0.1*rho_l*rho_0*(1./(2.*pi*1./50))**0.25*np.cos(2.*pi*(rr_lag-0.25)+2.*pi*(ss_lag-0.75))*np.exp(0.25/(1./50)*-((rr_lag-0.25)**2+(ss_lag-0.75)**2))
    myrho_0 = 1.0
    #myrho_0 = 10.0
    mymu_0 = 0.01
    mymu_0 = 0.0
    #Rho = myrho_0*np.exp(-2*pi*(ss_lag-0.5)/C_h)+lagregion*myrho_0*np.exp(-2*pi*(ss_lag-0.5)/C_h)*(1j/(wave_omega*C_h)*C_v*np.exp(2*pi*(ss_lag-0.5)/(2*C_h))*np.exp(1j*2*pi*(wave_k*(rr_lag-0.5)+wave_l*(ss_lag-0.5)))).real
    Rho = myrho_0*np.exp(-(ss_lag-0.5)/C_h)+lagregion*myrho_0*np.exp(-(ss_lag-0.5)/C_h)*(1j/(wave_omega*C_h)*C_v*np.exp((ss_lag-0.5)/(2*C_h))*np.exp(1j*(wave_k*(rr_lag-0.5)+wave_l*(ss_lag-0.5)))).real
    #Rho = myrho_0*np.exp(-2*pi*(ss_lag-0.5)/C_h)
    #Mu = mu_l*(mymu_0*np.exp(-2*pi*(ss_lag-0.5)/C_h)+lagregion*mymu_0*(1j/(wave_omega*C_h)*np.exp(-2*pi*(ss_lag-0.5)/C_h)*C_v*np.exp(2*pi*(ss_lag-0.5)/(2*C_h))*np.exp(1j*2*pi*(wave_k*(rr_lag-0.5)+wave_l*(ss_lag-0.5)))).real)
    Mu = mu_l*(mymu_0*np.exp(-(ss_lag-0.5)/C_h)+lagregion*mymu_0*(1j/(wave_omega*C_h)*np.exp(-(ss_lag-0.5)/C_h)*C_v*np.exp((ss_lag-0.5)/(2*C_h))*np.exp(1j*(wave_k*(rr_lag-0.5)+wave_l*(ss_lag-0.5)))).real)
    Mu = Mu.reshape(b_q,1).copy()
    ##make viscosity big outside of a circle of radius 0.5
    #Mu[np.nonzero(X_1**2+X_2**2 <= 0.25)[0]] = 0.01
    Rho = Rho.reshape(b_q,1).copy()
    rho_0 = np.min(Rho)
    #print rho_0
    Rho = Rho-np.min(Rho)
    
    ##surface waves
    #C_v = 0.01
    #C_v = 0.001
    #C_v = 0.0002
    #C_v = 0.00001
    #C_v = 0.00002
    ##C_v = 0.00001
    #C_v = 0.00004
    
    ##scale down C_v
    #scale_fac = 32
    #scale_fac = 1
    #C_v = C_v/scale_fac

    ##C_v = 0.0000025
    #dens_1 = 0.9
    #dens_0 = 1.1
    ##dens_1 = 0.999
    ##dens_0 = 1.001
    #rho_0 = dens_1
    #wave_k = 2.0
    #wave_omega = np.sqrt(grav*wave_k*(dens_0-dens_1)/(dens_0+dens_1))
    #wave_omega = np.sqrt(grav/(2*pi)*wave_k*np.tanh(pi*wave_k)*(dens_0-dens_1)/(dens_0+dens_1))
    #print 'period is', 1.0/wave_omega 
    #region0 = np.zeros((N*N,1))
    #xx_2_c_flat = xx_2_c.reshape(N*N,1)
    #xx_1_stag_flat = xx_1_stag.reshape(N*N,1)
    #xx_2_stag_flat = xx_2_stag.reshape(N*N,1)
    #xx_1_c_flat = xx_1_c.reshape(N*N,1)
    #region0[xx_2_c_flat < C_v*(np.sinh(pi*wave_k)*wave_k/(2*pi*1j*wave_omega)*np.exp(1j*(2*pi*wave_k*xx_1_stag_flat))).real] = 1.0
    #region0 = region0.flatten(1)
    #lagregion0 = np.zeros((b_q,1))
    #lagregion0[X_2 < C_v*(np.sinh(pi*wave_k)*wave_k/(2*pi*1j*wave_omega)*np.exp(1j*(2*pi*wave_k*(rr_lag.reshape(b_q,1)-0.5)))).real] = 1.0
    #lagregion0 = lagregion0.flatten(1)
    #Rho = np.zeros((b_q,1)) 
    #Rho[lagregion0 == 1] = dens_0-dens_1
    #Rho = Rho.reshape(b_q,1).copy()
    #Mu = np.zeros((b_q,1)) 
    #Mu = Mu.reshape(b_q,1).copy()
    #u_1_soln = np.zeros((N*N,1))
    #u_2_soln = np.zeros((N*N,1))
    #u_1_soln[region0 == 1] = C_v*(1j*wave_k*np.cosh(2*pi*wave_k*(xx_2_c_flat[region0 == 1]+0.5))*np.exp(1j*(wave_k*2*pi*xx_1_stag_flat[region0 == 1]))).real
    #u_2_soln[region0 == 1] = C_v*(wave_k*np.sinh(2*pi*wave_k*(xx_2_c_flat[region0 == 1]+0.5))*np.exp(1j*(wave_k*2*pi*xx_1_stag_flat[region0 == 1]))).real
    #u_1_soln[region0 == 0] = -C_v*(1j*wave_k*np.cosh(-2*pi*wave_k*(xx_2_stag_flat[region0 == 0]-0.5))*np.exp(1j*(wave_k*2*pi*xx_1_c_flat[region0 == 0]))).real
    #u_2_soln[region0 == 0] = C_v*(wave_k*np.sinh(-2*pi*wave_k*(xx_2_stag_flat[region0 == 0]-0.5))*np.exp(1j*(wave_k*2*pi*xx_1_c_flat[region0 == 0]))).real
    #u_1_soln = u_1_soln.reshape(N,N)
    #u_2_soln = u_2_soln.reshape(N,N)
    #u_1_old_soln = np.zeros((N*N,1))
    #u_2_old_soln = np.zeros((N*N,1))
    #u_1_old_soln[region0 == 1] = C_v*(1j*wave_k*np.cosh(2*pi*wave_k*(xx_2_c_flat[region0 == 1]+0.5))*np.exp(2*pi*1j*(wave_k*xx_1_stag_flat[region0 == 1]-wave_omega*dt))).real
    #u_2_old_soln[region0 == 1] = C_v*(wave_k*np.sinh(2*pi*wave_k*(xx_2_c_flat[region0 == 1]+0.5))*np.exp(2*pi*1j*(wave_k*xx_1_stag_flat[region0 == 1]-wave_omega*dt))).real
    #u_1_old_soln[region0 == 0] = -C_v*(1j*wave_k*np.cosh(-2*pi*wave_k*(xx_2_stag_flat[region0 == 0]-0.5))*np.exp(2*pi*1j*(wave_k*xx_1_c_flat[region0 == 0]-wave_omega*dt))).real
    #u_2_old_soln[region0 == 0] = C_v*(wave_k*np.sinh(-2*pi*wave_k*(xx_2_stag_flat[region0 == 0]-0.5))*np.exp(2*pi*1j*(wave_k*xx_1_c_flat[region0 == 0]-wave_omega*dt))).real
    #u_1_old_soln = u_1_old_soln.reshape(N,N)
    #u_2_old_soln = u_2_old_soln.reshape(N,N)

    ##make triangulated region circular
    #N_circle = 3*N_s
    #dcircle = 1./N_circle
    #s_circle = np.arange(0+dcircle/2, 1.0+dcircle/2, dcircle)
    #X_cen = -0.5+myratio*0.5 
    #Y_cen = -0.5+myratio*0.5 
    #X_circle = X_cen+0.5*myratio*np.sqrt(2)*np.cos(2*pi*s_circle)
    #Y_circle = Y_cen+0.5*myratio*np.sqrt(2)*np.sin(2*pi*s_circle)
    #X_circle = X_circle.reshape(N_circle,1)
    #Y_circle = Y_circle.reshape(N_circle,1)
    #r_out_lag = np.arange(dr/2, np.sqrt(2)*1.0+dr/2, dr)
    #s_out_lag = np.arange(ds/2, np.sqrt(2)*1.0+ds/2, ds)
    #ss_out_lag, rr_out_lag = np.meshgrid(s_out_lag, r_out_lag)
    #N_out = len(ss_out_lag.flatten(1))
    ##rr_out_lag = rr_out_lag.flatten(1)
    ##ss_out_lag = ss_out_lag.flatten(1)
    #rr_out_lag = rr_out_lag.reshape(N_out,1).flatten(1)
    #ss_out_lag = ss_out_lag.reshape(N_out,1).flatten(1)

    #X_outside = (myratio*rr_out_lag-0.5-0.5*myratio*(np.sqrt(2)-1)).copy()
    #Y_outside = (myratio*ss_out_lag-0.5-0.5*myratio*(np.sqrt(2)-1)).copy()
    #X_outside_temp = X_outside[(((X_outside-X_cen)**2+(Y_outside-Y_cen)**2) < 0.5*myratio**2) * ~((X_outside > -0.5) * (X_outside < (myratio-0.5)) * (Y_outside >  -0.5) * (Y_outside < (myratio-0.5)))][:,np.newaxis]
    #Y_outside = Y_outside[(((X_outside-X_cen)**2+(Y_outside-Y_cen)**2) < 0.5*myratio**2) * ~((X_outside > -0.5) * (X_outside < (myratio-0.5)) * (Y_outside >  -0.5) * (Y_outside < (myratio-0.5)))][:,np.newaxis]
    #X_outside = X_outside_temp.copy()
    #X_1 = np.vstack((X_1,X_circle,X_outside))
    #X_2 = np.vstack((X_2,Y_circle,Y_outside))
    ##X_1 = X_outside
    ##X_2 = Y_outside

    #edges = round_trip_connect(b_q, b_q+N_circle-1) 
    #edges_array = np.asarray(edges)
    #Mu = np.vstack((Mu,np.zeros((X_1.shape[0]-b_q,1))))
    #Rho = np.vstack((Rho,np.zeros((X_1.shape[0]-b_q,1))))
    #b_q = X_1.shape[0] 
    #print 'number of Lagrangian points is', b_q
    #info = triangle.MeshInfo()
    #X_arr = np.hstack((X_1,X_2)).copy()
    #info.set_points(list(X_arr), range(2,b_q+2))
    #info.set_facets(edges)
    #mesh = triangle.build(info, allow_volume_steiner=False,
    #   quality_meshing=False, allow_boundary_steiner=False)
    #N_tri = len(mesh.elements) 
    #TV = np.array(mesh.elements).copy()
    #markers = np.array(mesh.point_markers)
    #inds = np.nonzero(markers <= 1)[0] 
    #if (len(inds) != 0):
    #    print 'added new point illegally'
    #    sys.exit()

    if (np.mod(N_tri, 16) == 0):
        b_WGPSIZE_3 = 16
    elif (np.mod(N_tri, 8) == 0):
        b_WGPSIZE_3 = 8
    elif (np.mod(N_tri, 4) == 0):
        b_WGPSIZE_3 = 4
    elif (np.mod(N_tri, 2) == 0):
        b_WGPSIZE_3 = 2
    else:
        b_WGPSIZE_3 = 1
    print 'new number of triangles is', N_tri

    Rho_tot = 0.0
    for q in range(0,N_tri):
        verts = TV[q,:]
        X11 = X_1[verts[0]]
        X12 = X_2[verts[0]]
        X21 = X_1[verts[1]]
        X22 = X_2[verts[1]]
        X31 = X_1[verts[2]]
        X32 = X_2[verts[2]]

        a1 = X21-X11
        a2 = X22-X12
        b1 = X31-X11
        b2 = X32-X12

        norm = a1*b2-a2*b1
        if (norm < 0):
            print 'negative normal'
            print 'im here'
            temp = TV[q,0]
            TV[q,0] = TV[q,1]
            TV[q,1] = temp
        Rho_avg = 1.0/3.0*(Rho[verts[0]]+Rho[verts[1]]+Rho[verts[2]])
        Rho_tot = Rho_tot + 0.5*np.abs(norm)*Rho_avg 

    rho_av=(Rho_tot+dx**2*rho_0*N**2)/(dx**2*N**2);
    print 'rho_av is', rho_av

    mu_min = mu_0
    mu_max = mu_0 + mu_0*mu_l
    #mu_hat = 0.5*(mu_min+mu_max)
    rho_min = rho_0
    rho_max = rho_0 + rho_0*rho_l
    rho_max = rho_0+np.max(Rho)
    mu_max = mu_0+np.max(Mu)
    #rho_hat = 0.5*(rho_min+rho_max)
    mu_hat = mu_max
    rho_hat = rho_max

    #rho_hat = rho_0

    #mu_bar = dtheta*np.max(Mu)
    #rho_bar = dtheta*np.max(Rho)
    #mu_bar = 0.025 
    #rho_bar = 0.5 
    #mu_bar = 10.0*mu_0 
    #rho_bar = 10.0*rho_0 
    #mu_bar = mu_0
    #rho_bar = 0.5*rho_0
    #mu_bar = 2*mu_0
    #rho_bar = 100*rho_0
    #mu_bar = 0.0
    #rho_bar = 0.0
    #print 'mu_bar is', mu_bar
    #print 'rho_bar is', rho_bar
    #mu_0 = mu_0+mu_bar
    #rho_0 = rho_0+rho_bar

    
    ##include for tethering
    N_teth = 2*N
    #N_teth = 4*N-1
    #N_teth = N
    dx_teth = dx/2.0
    #dx_teth = dx
    #X_1_teth = np.arange(-0.5+dx_teth/2, 0.5+dx_teth/2, dx_teth)
    X_1_teth = np.arange(-0.5, 0.5, dx_teth)
    #X_1_teth = X_1_teth.reshape(N_teth,1).copy()
    X_1_teth = X_1_teth.reshape(2*N,1).copy()
    #X_1_teth = np.vstack((X_1_teth,(-0.5)*np.ones((2*N-1,1))))
    #X_2_teth = (-0.5+dx/2)*np.ones((N_teth,1))
    #X_2_teth = (-0.5)*np.ones((N_teth,1))
    X_2_teth = (-0.5)*np.ones((2*N,1))
    #X_2_teth = np.vstack((X_2_teth,(np.arange(-0.5+dx_teth/2, 0.5+dx_teth/2, dx_teth)).reshape(N_teth/2,1)))
    #X_2_teth = np.vstack((X_2_teth,(np.arange(-0.5+dx_teth, 0.5, dx_teth)).reshape(2*N-1,1)))
    k_teth = 1.0*10**-12 
    k_teth = 1.0*10**3
    k_teth = 1.0*10**4
    #k_teth = 1.0
    #k_teth = 1.0*10**7
    #k_teth = 4.0*10**5
    #k_teth = 4.0*10**4
    #k_teth = 1.0*10**6
    #k_teth = 0.0
    #k_teth = 1.0*10**1
    #k_teth = 4.0*10**3
    #k_teth = 4.0*10**8
    #k_teth = 1.0*10**1
    #N_teth = 2*N_teth
    #X_1_teth_old = X_1_teth.copy()
    #X_1_teth = np.vstack((X_1_teth, X_2_teth))
    #X_2_teth = np.vstack((X_2_teth, X_1_teth_old))
    X_1 = np.vstack((X_1, X_1_teth))
    X_2 = np.vstack((X_2, X_2_teth))
    X_2_save = X_2
    if (np.mod(N_teth, 8) == 0):
        b_WGPSIZE_4 = 8
    elif (np.mod(N_teth, 4) == 0):
        b_WGPSIZE_4 = 4
    elif (np.mod(N_teth, 2) == 0):
        b_WGPSIZE_4 = 2
    else:
        b_WGPSIZE_4 = 1
    if (np.mod(b_q+N_teth, 16) == 0):
        b_WGPSIZE_2 = 16
    if (np.mod(b_q+N_teth, 8) == 0):
        b_WGPSIZE_2 = 8
    elif (np.mod(b_q+N_teth, 4) == 0):
        b_WGPSIZE_2 = 4
    elif (np.mod(b_q+N_teth, 2) == 0):
        b_WGPSIZE_2 = 2
    else:
        b_WGPSIZE_2 = 1

    firstn1 = -1*np.ones((N,N))
    nextn1 = -1*np.ones((b_q+N_teth,1))
    firstn2 = -1*np.ones((N,N))
    nextn2 = -1*np.ones((b_q+N_teth,1))
    numn1 = np.zeros((N,N))
    numn2 = np.zeros((N,N))
    N_pre_teth = b_q
    #b_q = b_q+N_teth

    Mu = np.vstack((Mu, np.zeros((N_teth,1))))
    Rho = np.vstack((Rho, np.zeros((N_teth,1))))

    start_time = time.time()
    for k in range(0,N_teth):
        s_1_stag = (X_1[N_pre_teth+k]+0.5)/dx
        s_2_stag = (X_2[N_pre_teth+k]+0.5)/dx        
        i_1_stag = (np.floor(s_1_stag)).astype(np.uint32)
        i_2_stag = (np.floor(s_2_stag)).astype(np.uint32)
        s_1_c = (X_1[N_pre_teth+k]+0.5-dx/2)/dx;
        s_2_c = (X_2[N_pre_teth+k]+0.5-dx/2)/dx;        
        i_1_c = (np.remainder(np.floor(s_1_c),N)).astype(np.uint32);
        i_2_c = (np.remainder(np.floor(s_2_c),N)).astype(np.uint32);
        if (firstn1[i_2_c,i_1_stag]==-1):
            firstn1[i_2_c,i_1_stag]=N_pre_teth+k 
            numn1[i_2_c,i_1_stag] += 1
        else:
            nextn1[N_pre_teth+k]=firstn1[i_2_c,i_1_stag]
            firstn1[i_2_c,i_1_stag]=N_pre_teth+k
            numn1[i_2_c,i_1_stag] += 1
        if (firstn2[i_2_stag,i_1_c]==-1):
            firstn2[i_2_stag,i_1_c]=N_pre_teth+k 
            numn2[i_2_stag,i_1_c] += 1
        else:
            nextn2[N_pre_teth+k]=firstn2[i_2_stag,i_1_c]
            firstn2[i_2_stag,i_1_c]=N_pre_teth+k
            numn2[i_2_stag,i_1_c] += 1
    print 'time spent binning points is', time.time()-start_time, 'seconds'

    print 'numn1 is', np.sum(numn1)
    print 'numn2 is', np.sum(numn2)

    L_tf = calc_transf(dx)
    i_1, i_2, ip_1, ip_2, im_1, im_2 = calc_inds(N)

    #u_1_soln = 2*np.cos(2*pi*xx_1_stag)*np.sin(4*pi*xx_2_c)
    #u_2_soln = -np.sin(2*pi*xx_1_c)*np.cos(4*pi*xx_2_stag)
    #u_1_soln = np.zeros((N,N)) #initial conditions
    #u_2_soln = np.zeros((N,N))
    #u_2_soln = -0.1*rho_0*rho_l*np.exp(-0.5*(xx_2_stag+0.5))*np.sin(4*pi*(-xx_1_c-0.5)-pi/2)*np.sin(2*pi*(-xx_2_stag-0.5)-pi/2) 
    #u_1_soln = np.sin(2*pi*xx_2_c) #shear field
    #u_2_soln = np.zeros((N,N))

    #move data to gpu
    start_time = time.time()
    plan = Plan((N, N), queue=queue, dtype=np.complex128, fast_math=False)
    print 'time spent making plan is', time.time()-start_time, 'seconds'

    u_1_gpu = cl_array.to_device(queue, u_1_soln.astype(np.float64))
    u_2_gpu = cl_array.to_device(queue, u_2_soln.astype(np.float64))
    #u_1_o_gpu = cl_array.to_device(queue, u_1_soln.astype(np.float64))
    #u_2_o_gpu = cl_array.to_device(queue, u_2_soln.astype(np.float64))
    u_1_o_gpu = cl_array.to_device(queue, u_1_old_soln.astype(np.float64))
    u_2_o_gpu = cl_array.to_device(queue, u_2_old_soln.astype(np.float64))
    u_1_tf_gpu = cl_array.to_device(queue, u_1_soln.astype(np.complex128))
    u_2_tf_gpu = cl_array.to_device(queue, u_2_soln.astype(np.complex128))
    L_tf_gpu = cl_array.to_device(queue, L_tf.astype(np.complex128))
    firstn1_gpu = cl_array.to_device(queue, firstn1.astype(np.int32))
    nextn1_gpu = cl_array.to_device(queue, nextn1.astype(np.int32))
    firstn2_gpu = cl_array.to_device(queue, firstn2.astype(np.int32))
    nextn2_gpu = cl_array.to_device(queue, nextn2.astype(np.int32))
    f_1_gpu = cl_array.empty_like(u_1_gpu) 
    f_2_gpu = cl_array.empty_like(u_1_gpu) 
    f_1_tf_gpu = cl_array.empty_like(u_1_tf_gpu) 
    f_2_tf_gpu = cl_array.empty_like(u_1_tf_gpu) 
    mu_1_gpu = cl_array.empty_like(u_1_gpu) 
    mu_2_gpu = cl_array.empty_like(u_1_gpu) 
    if SAVE:
        rho_1_save_gpu = cl_array.empty_like(u_1_gpu) 
        rho_2_save_gpu = cl_array.empty_like(u_1_gpu) 
        mu_1_save_gpu = cl_array.empty_like(u_1_gpu) 
        mu_2_save_gpu = cl_array.empty_like(u_1_gpu) 
    if PR_INT:
        vort_gpu = cl_array.empty_like(u_1_gpu) 
    rho_1_gpu = cl_array.empty_like(u_1_gpu) 
    rho_2_gpu = cl_array.empty_like(u_1_gpu) 
    p_gpu = cl_array.empty_like(u_1_gpu) 
    p_last_gpu = cl_array.empty_like(u_1_gpu) 
    p_tf_gpu = cl_array.empty_like(u_1_tf_gpu) 
    u_1_tilde_gpu = cl_array.empty_like(u_1_gpu) 
    u_2_tilde_gpu = cl_array.empty_like(u_1_gpu) 
    X_1_gpu = cl_array.to_device(queue, X_1.astype(np.float64))
    X_2_gpu = cl_array.to_device(queue, X_2.astype(np.float64))
    X_2_save_gpu = cl_array.to_device(queue, X_2.astype(np.float64))
    X_1_tilde_gpu = cl_array.to_device(queue, X_1.astype(np.float64))
    X_2_tilde_gpu = cl_array.to_device(queue, X_2.astype(np.float64))
    rr_lag_gpu = cl_array.to_device(queue, rr_lag.astype(np.float64))
    F_1_gpu = cl_array.empty_like(X_1_gpu) 
    F_2_gpu = cl_array.empty_like(X_1_gpu) 
    Mu_gpu = cl_array.to_device(queue, Mu.astype(np.float64)) 
    Rho_gpu = cl_array.to_device(queue, Rho.astype(np.float64)) 
    U_1_gpu = cl_array.empty_like(X_1_gpu) 
    U_2_gpu = cl_array.empty_like(X_1_gpu) 
    U_2_gpu = cl_array.empty_like(X_1_gpu) 
    region_gpu = cl_array.to_device(queue, region.astype(np.float64))

    #for triangulation
    TV_1_gpu = cl_array.to_device(queue, TV[:,0].astype(np.uint32))
    TV_2_gpu = cl_array.to_device(queue, TV[:,1].astype(np.uint32))
    TV_3_gpu = cl_array.to_device(queue, TV[:,2].astype(np.uint32))
    en = np.zeros((1,1))
    en_gpu = cl_array.to_device(queue, en.astype(np.float64))
    en_sum = np.zeros((N_tri,1))
    en_sum_gpu = cl_array.to_device(queue, en_sum.astype(np.float64))

    #for tethering
    X_1_teth_gpu = cl_array.to_device(queue, X_1_teth.astype(np.float64))
    X_2_teth_gpu = cl_array.to_device(queue, X_2_teth.astype(np.float64))

    #ad hoc - put in d'Alembert form eventually
    rho_1_extra_gpu = cl_array.empty_like(u_1_gpu) 
    rho_2_extra_gpu = cl_array.empty_like(u_1_gpu) 

    ##add forcing term
    #xx_1_stag_gpu = cl_array.to_device(ctx, queue, xx_1_stag.astype(np.float64))
    #xx_2_stag_gpu = cl_array.to_device(ctx, queue, xx_2_stag.astype(np.float64))
    #xx_1_c_gpu = cl_array.to_device(ctx, queue, xx_1_c.astype(np.float64))
    #xx_2_c_gpu = cl_array.to_device(ctx, queue, xx_2_c.astype(np.float64))

    #solve
    queue.finish()
    start_time = time.time()
    #retriang(0)

    for t_step in range(1, clock_max+1):
        solve(t_step)

        if (t_step-1) % DISP_INT == 0:
            queue.finish()
            print t_step
            print 'time is', t_step*dt
            if SAVE: 
                #print u_1_gpu.get().max()
                print 'u_1 max is', u_1_gpu.get().max()
                print 'u_2 max is', u_2_gpu.get().max()
                print 'rho_1 max is', rho_1_save_gpu.get().max()
                print 'rho_1 min is', rho_1_save_gpu.get().min()
                print 'mu_1 max is', mu_1_save_gpu.get().max()
                print 'mu_1 min is', mu_1_save_gpu.get().min()
                print 'av height is', np.mean(X_2_gpu.get()[N_r:b_q])
                print 'mom1 is', (u_1_gpu.get()*(rho_1_save_gpu.get())).sum()*dx*dx
                print 'mom2 is', (u_2_gpu.get()*(rho_2_save_gpu.get())).sum()*dx*dx
                #f.write('%g \n' % np.mean(X_2_gpu.get()[0:N_r]))
        if PR_INT:
            if (t_step-1) % PR_INT == 0:
                #X_1 = X_1_gpu.get()
                #X_2 = X_2_gpu.get()
                #for q in range(0,N_tri):
                #    verts = TV[q,:]
                #    X11 = X_1[verts[0]]
                #    X12 = X_2[verts[0]]
                #    X21 = X_1[verts[1]]
                #    X22 = X_2[verts[1]]
                #    X31 = X_1[verts[2]]
                #    X32 = X_2[verts[2]]

                #    a1 = X21-X11
                #    a2 = X22-X12
                #    b1 = X31-X11
                #    b2 = X32-X12

                #    norm = a1*b2-a2*b1
                #    if (norm < 0):
                #        print 'negative normal'

                from os.path import join as pjoin
                #vis_p = p_gpu.get().reshape(N,N,1)
                #file_name = "p-%04d.vts"%t_step
                #path_to_file = pjoin("vis_files", file_name) 
                #write_structured_grid(path_to_file, mesh3.reshape(3, N, N, 1), 
                #point_data=[("p", vis_p)])

                #file_name = "pslice-%04d.txt"%t_step
                #vis_p = p_gpu.get()
                #vis_p = vis_p[N/2,:]
                #path_to_file = pjoin("vis_files", file_name) 
                #np.savetxt(path_to_file, vis_p)

                #vis_mu = mu_0+mu_1_save_gpu.get().reshape(N,N,1)
                vis_u_1 = u_1_gpu.get().reshape(N,N,1)
                file_name = "u_1-%04d.vts"%t_step
                path_to_file = pjoin("vis_files", file_name) 
                #write_structured_grid(path_to_file, mesh3.reshape(3, N, N, 1), 
                #point_data=[("u_1", vis_u_1)])

                #vis_u_2 = region*myrho_0*np.exp(-xx_2_stag/C_h)*(1j/(wave_omega*C_h)*C_v*np.exp(xx_2_stag/(2*C_h))*np.exp(1j*(wave_k*xx_1_c+wave_l*xx_2_stag))).real
                vis_u_2 = u_2_gpu.get().reshape(N,N,1)
                file_name = "u_2-%04d.vts"%t_step
                path_to_file = pjoin("vis_files", file_name) 
                #write_structured_grid(path_to_file, mesh3.reshape(3, N, N, 1), 
                #point_data=[("u_2", vis_u_2)])

                #vis_u_2 = u_2_gpu.get().reshape(N,N,1)
                #file_name = "u_2-%04d.vts"%t_step
                #path_to_file = pjoin("vis_files", file_name) 
                #write_structured_grid(path_to_file, mesh3.reshape(3, N, N, 1), 
                #point_data=[("u_2", vis_u_2)])

                #vis_mu_2 = mu_2_save_gpu.get().reshape(N,N,1)
                #vis_mu_2 = (rho_2_save_gpu.get()-1.0).reshape(N,N,1)
                #if t_step == 1:
		#    myavm = np.mean(sqregion*(rho_2_save_gpu.get()-(myrho_0*np.exp(-2*pi*xx_2_stag/C_h))))
		#vis_mu_2 = ((1-sqregion)*myavm+sqregion*(rho_2_save_gpu.get()-(myrho_0*np.exp(-2*pi*xx_2_stag/C_h)))).reshape(N,N,1) #subtract background
		#vis_mu_2 = ((rho_2_save_gpu.get()-(myrho_0*np.exp(-2*pi*xx_2_stag/C_h)))).reshape(N,N,1) #subtract background
		vis_mu_2 = ((rho_2_save_gpu.get()-(myrho_0*np.exp(-xx_2_stag/C_h)))).reshape(N,N,1) #subtract background

                ##scale up mu_2
                #vis_mu_2 = np.zeros((N,N))
                #for i in range(0,(N/2)/scale_fac):
                #    for j in range(0,scale_fac):
                #        vis_mu_2[N/2-1+scale_fac*i+j,:] = ((1.0-j/scale_fac)*rho_2_save_gpu.get()[N/2-1+i,:]+j/scale_fac*rho_2_save_gpu.get()[N/2+i,:])-1.0
                #        vis_mu_2[N/2-1-scale_fac*i-j,:] = ((1.0-j/scale_fac)*rho_2_save_gpu.get()[N/2-1-i,:]+j/scale_fac*rho_2_save_gpu.get()[N/2-2-i,:])-1.0
                #vis_mu_2 = vis_mu_2.reshape(N,N,1)

		#myav = np.mean(sqregion*(rho_2_save_gpu.get()-(myrho_0*np.exp(-2*pi*xx_2_stag/C_h)))) #subtract background
		#vis_mu_2 = ((1-sqregion)*myav+sqregion*(rho_2_save_gpu.get()-(myrho_0*np.exp(-2*pi*xx_2_stag/C_h)))).reshape(N,N,1) #subtract background
		#vis_mu_2 = ((rho_2_save_gpu.get()-rho_0-(myrho_0*np.exp(-2*pi*xx_2_stag/C_h)))).reshape(N,N,1) #subtract background
		#vis_mu_2 = (sqregion*(mu_2_save_gpu.get()-mu_0-mu_l*(mymu_0*np.exp(-2*pi*xx_2_stag/C_h)))).reshape(N,N,1) #subtract background
                file_name = "mu_2-%04d.vts"%t_step
                path_to_file = pjoin("vis_files", file_name) 
                write_structured_grid(path_to_file, mesh3.reshape(3, N, N, 1), 
                point_data=[("mu_2", vis_mu_2)])

                vis_rho = rho_1_save_gpu.get().reshape(N,N,1)
                #vis_en = (rho_1_save_gpu.get()*(u_1_gpu.get()*u_1_gpu.get()+u_2_gpu.get()*u_2_gpu.get())).reshape(N,N,1)

                #if t_step == 1:
                #    myavs = np.mean(sqregion*0.5*((myrho_0*np.exp(-2*pi*xx_2_c/C_h))*u_1_gpu.get()*u_1_gpu.get()+(myrho_0*np.exp(-2*pi*xx_2_stag/C_h))*u_2_gpu.get()*u_2_gpu.get())) 

                #vis_en = (myavs*(1-sqregion)+sqregion*0.5*((myrho_0*np.exp(-2*pi*xx_2_c/C_h))*u_1_gpu.get()*u_1_gpu.get()+(myrho_0*np.exp(-2*pi*xx_2_stag/C_h))*u_2_gpu.get()*u_2_gpu.get())).reshape(N,N,1)
                #vis_en = (0.5*((myrho_0*np.exp(-2*pi*xx_2_c/C_h))*u_1_gpu.get()*u_1_gpu.get()+(myrho_0*np.exp(-2*pi*xx_2_stag/C_h))*u_2_gpu.get()*u_2_gpu.get())).reshape(N,N,1)
                vis_en = (0.5*((myrho_0*np.exp(-xx_2_c/C_h))*u_1_gpu.get()*u_1_gpu.get()+(myrho_0*np.exp(-xx_2_stag/C_h))*u_2_gpu.get()*u_2_gpu.get())).reshape(N,N,1)
                #vis_en = (0.5*((myrho_0*np.exp(-xx_2_c/C_h))*(1j/wave_k*(1j*wave_l+1.0/(2*C_h))*C_v*np.exp(xx_2_c/(2*C_h))*np.exp(1j*(wave_k*xx_1_stag+wave_l*xx_2_c))).real*region*(1j/wave_k*(1j*wave_l+1.0/(2*C_h))*C_v*np.exp(xx_2_c/(2*C_h))*np.exp(1j*(wave_k*xx_1_stag+wave_l*xx_2_c))).real*region+(myrho_0*np.exp(-xx_2_stag/C_h))*C_v*np.exp(xx_2_stag/(2*C_h))*(np.exp(1j*(wave_k*xx_1_c+wave_l*xx_2_stag))).real*region*C_v*np.exp(xx_2_stag/(2*C_h))*(np.exp(1j*(wave_k*xx_1_c+wave_l*xx_2_stag))).real*region)).reshape(N,N,1)
                print 'kinetic energy is', np.sum(vis_en)*dx*dx
                #u_2_soln = C_v*np.exp(xx_2_stag/(2*C_h))*(np.exp(1j*(wave_k*xx_1_c+wave_l*xx_2_stag))).real*region
                #u_1_soln = (1j/wave_k*(1j*wave_l+1.0/(2*C_h))*C_v*np.exp(xx_2_c/(2*C_h))*np.exp(1j*(wave_k*xx_1_stag+wave_l*xx_2_c))).real*region

                #myav = np.mean(sqregion*0.5*grav*grav/(myrho_0*np.exp(-2*pi*xx_2_c/C_h)*buoy_freq*buoy_freq)*(rho_1_save_gpu.get()-(myrho_0*np.exp(-2*pi*xx_2_c/C_h)))*(rho_1_save_gpu.get()-(myrho_0*np.exp(-2*pi*xx_2_c/C_h))))

                #if t_step == 1:
                #    myavr = np.mean(sqregion*0.5*grav*C_h/(2*pi)/(myrho_0*np.exp(-2*pi*xx_2_c/C_h))*(rho_1_save_gpu.get()-(myrho_0*np.exp(-2*pi*xx_2_c/C_h)))*(rho_1_save_gpu.get()-(myrho_0*np.exp(-2*pi*xx_2_c/C_h))))
                #add potential energy
                #vis_en = vis_en + (myav*(1-sqregion)+sqregion*(0.5*grav*grav/(myrho_0*np.exp(-2*pi*xx_2_c/C_h)*buoy_freq*buoy_freq)*(rho_1_save_gpu.get()-(myrho_0*np.exp(-2*pi*xx_2_c/C_h)))*(rho_1_save_gpu.get()-(myrho_0*np.exp(-2*pi*xx_2_c/C_h))))).reshape(N,N,1)
                #vis_en = vis_en + (myavr*(1-sqregion)+sqregion*0.5*grav*C_h/(2*pi)/(myrho_0*np.exp(-2*pi*xx_2_c/C_h))*(rho_1_save_gpu.get()-(myrho_0*np.exp(-2*pi*xx_2_c/C_h)))*(rho_1_save_gpu.get()-(myrho_0*np.exp(-2*pi*xx_2_c/C_h)))).reshape(N,N,1)
                #vis_en = vis_en + (0.5*grav*C_h/(2*pi)/(myrho_0*np.exp(-2*pi*xx_2_c/C_h))*(rho_1_save_gpu.get()-(myrho_0*np.exp(-2*pi*xx_2_c/C_h)))*(rho_1_save_gpu.get()-(myrho_0*np.exp(-2*pi*xx_2_c/C_h)))).reshape(N,N,1)
                #vis_en2 = (0.5*grav*C_h/(2*pi)/(myrho_0*np.exp(-2*pi*xx_2_c/C_h))*(rho_1_save_gpu.get()-(myrho_0*np.exp(-2*pi*xx_2_c/C_h)))*(rho_1_save_gpu.get()-(myrho_0*np.exp(-2*pi*xx_2_c/C_h)))).reshape(N,N,1)
                vis_en2 = (0.5*grav*C_h/(myrho_0*np.exp(-xx_2_c/C_h))*(rho_1_save_gpu.get()-(myrho_0*np.exp(-xx_2_c/C_h)))*(rho_1_save_gpu.get()-(myrho_0*np.exp(-xx_2_c/C_h)))).reshape(N,N,1)
                #vis_en = vis_en + (0.5*grav*C_h/(myrho_0*np.exp(-xx_2_c/C_h))*myrho_0*region*np.exp(-xx_2_c/C_h)*(1j/(wave_omega*C_h)*C_v*np.exp(xx_2_c/(2*C_h))*np.exp(1j*(wave_k*xx_1_stag+wave_l*xx_2_c))).real*region*myrho_0*np.exp(-xx_2_c/C_h)*(1j/(wave_omega*C_h)*C_v*np.exp(xx_2_c/(2*C_h))*np.exp(1j*(wave_k*xx_1_stag+wave_l*xx_2_c))).real).reshape(N,N,1)
                print 'potential energy is', np.sum(vis_en2)*dx*dx
                vis_en = vis_en+vis_en2
                print 'total energy is', np.sum(vis_en)*dx*dx

                #vis_en = (rho_1_save_gpu.get()-(myrho_0*np.exp(-2*pi*(xx_2_c)/C_h))).reshape(N,N,1)
                #vis_en = rho_1_save_gpu.get().reshape(N,N,1)

                #myav = np.mean(sqregion*grav*xx_2_c*(rho_1_save_gpu.get()-(myrho_0*np.exp(-2*pi*xx_2_c/C_h))))
                #vis_en = vis_en + (myav*(1-sqregion)+sqregion*grav*xx_2_c*(rho_1_save_gpu.get()-(myrho_0*np.exp(-2*pi*xx_2_c/C_h)))).reshape(N,N,1)
                #add_gpu(X_2_gpu,X_2_save_gpu,-1.0,X_2_tilde_gpu)
                #calc_coefs_kernel(queue, (N_tri, 1), (b_WGPSIZE_3,1), X_1_gpu.data, X_2_gpu.data, TV_1_gpu.data, TV_2_gpu.data, TV_3_gpu.data, X_2_tilde_gpu.data, Mu_gpu.data, dx, N, rho_2_gpu.data, mu_2_gpu.data, mu_1_gpu.data, mu_2_gpu.data) 
                #vis_en = vis_en + (0.5*buoy_freq*buoy_freq*rho_1_gpu.get()*rho_2_gpu.get()*rho_2_gpu.get()).reshape(N,N,1)
                #vis_en = grav*xx_2_c.reshape(N,N,1)
                #vis_en = (rho_1_save_gpu.get()-rho_av).reshape(N,N,1)
                #vis_en = grav*xx_2_c*(rho_1_save_gpu.get()-rho_av).reshape(N,N,1)
                file_name = "rho-%04d.vts"%t_step
                path_to_file = pjoin("vis_files", file_name) 
                write_structured_grid(path_to_file, mesh3.reshape(3, N, N, 1), 
                point_data=[("rho", vis_en)])
                #f.write('%g %g \n' % (np.sum(xx_1_stag*vis_en.reshape(N,N))/np.sum(vis_en.reshape(N,N)),np.sum(xx_2_c*vis_en.reshape(N,N))/np.sum(vis_en.reshape(N,N))) )

                #vis_vort = vort_gpu.get().reshape(N,N,1)
                #file_name = "vort-%04d.vts"%t_step
                #path_to_file = pjoin("vis_files", file_name) 
                #write_structured_grid(path_to_file, mesh3.reshape(3, N, N, 1), 
                #point_data=[("vort", vis_vort)])

                #vis_f_1 = f_1_gpu.get().reshape(N,N,1)
                #file_name = "f_1-%04d.vts"%t_step
                #path_to_file = pjoin("vis_files", file_name) 
                #write_structured_grid(path_to_file, mesh3.reshape(3, N, N, 1), 
                #point_data=[("f_1", vis_f_1)])

                #vis_f_2 = f_2_gpu.get().reshape(N,N,1)
                #file_name = "f_2-%04d.vts"%t_step
                #path_to_file = pjoin("vis_files", file_name) 
                #write_structured_grid(path_to_file, mesh3.reshape(3, N, N, 1), 
                #point_data=[("f_2", vis_f_2)])

                lims = np.array((-0.5,0.5))
                lims = lims.reshape(2,1).copy() 
                #vis_X_1 = np.vstack((np.mod(X_1_gpu.get()+0.5, 1)-0.5,lims))
                #vis_X_2 = np.vstack((np.mod(X_2_gpu.get()+0.5, 1)-0.5,lims))
                #vis_X_1 = np.vstack((X_1_gpu.get(),lims))
                #vis_X_2 = np.vstack((X_2_gpu.get(),lims))
                vis_X_1 = np.vstack((X_1_gpu.get()[0:b_q],lims))
                vis_X_2 = np.vstack((X_2_gpu.get()[0:b_q],lims))
                #vis_X_1 = vis_X_1[0:N_r+2].copy()
                #vis_X_2 = vis_X_2[0:N_r+2].copy()
                #print vis_X_1[0:5].copy().shape
                #print vis_X_1.shape
                #vis_F_1 = np.vstack((F_1_gpu.get().reshape(b_q,1),np.zeros((2,1))))
                #vis_F_2 = np.vstack((F_2_gpu.get().reshape(b_q,1),np.zeros((2,1))))
                vis_F_1 = np.vstack((F_1_gpu.get()[0:b_q].reshape(b_q,1),np.zeros((2,1))))
                vis_F_2 = np.vstack((F_2_gpu.get()[0:b_q].reshape(b_q,1),np.zeros((2,1))))
                #vis_F_2 = F_2_gpu.get().reshape(b_q).copy()
                #vis_F_1 = vis_F_1[0:N_r+2].copy()
                #vis_F_2 = vis_F_2[0:N_r+2].copy()
                #print vis_X_1.shape
                #print X_1_gpu.get().shape
                vis_X_3 = np.zeros((vis_X_1.shape))
                vis_X = np.hstack((vis_X_1,vis_X_2,vis_X_3)).copy()
                #print F_1_gpu.get().reshape(b_q).copy().shape
                data = [
                #("F_1", F_1_gpu.get().reshape(b_q+2).copy()),
                #("F_2", F_2_gpu.get().reshape(b_q+2).copy()),
                #("F_1", vis_F_1.reshape(b_q+2).copy()),
                #("F_2", vis_F_2.reshape(b_q+2).copy()),
                #("F_1", vis_F_1.copy()),
                #("F_2", vis_F_2.copy()),
                ]

                file_name = "X-%04d.vtu"%t_step
                compressor = None

                grid = UnstructuredGrid(
                     (b_q+2, DataArray("points", vis_X, vector_format=VF_LIST_OF_VECTORS)),
                     cells=np.arange(b_q+2, dtype=np.uint32),
                     cell_types=np.asarray([VTK_VERTEX] * (b_q+2), dtype=np.uint8))
                for name, field in data:
                   grid.add_pointdata(DataArray(name, field,
                   vector_format=VF_LIST_OF_COMPONENTS))

                from os.path import exists
                path_to_file = pjoin("vis_files", file_name)
                if exists(path_to_file):
                   raise RuntimeError("output file '%s' already exists"
                   % path_to_file)

                #outf = open(path_to_file, "w")
                #AppendedDataXMLGenerator(compressor)(grid).write(outf)
                #outf.close()

                #vis_X_1 = X_1_gpu.get()-np.mean(X_1_gpu.get())
                #vis_X_2 = X_2_gpu.get()-np.mean(X_2_gpu.get())
                #vis_X_3 = np.zeros((vis_X_1.shape))
                #vis_X = np.hstack((vis_X_1,vis_X_2,vis_X_3)).copy()
                #lims_tri = np.array([[b_q,b_q,b_q],[b_q+1,b_q+1,b_q+1]])
                ##lims_tri = np.array([[b_q,b_q,b_q]])
                #vis_tri = np.vstack((TV, lims_tri))
                vis_tri = TV

                vis_X_1 = X_1_gpu.get()
                vis_X_2 = X_2_gpu.get()
                vis_X_3 = np.zeros((vis_X_1.shape))
                vis_X = np.hstack((vis_X_1,vis_X_2,vis_X_3)).copy()

                file_name = "tri-%04d.vtu"%t_step
                compressor = None
                cell_data = []
                point_data = []

                grid = UnstructuredGrid(
                    (b_q+N_teth+b_add,
                    DataArray("points", vis_X,
                    vector_format=VF_LIST_OF_VECTORS)),
                    cells=vis_tri.reshape(-1),
                    cell_types=np.asarray([VTK_TRIANGLE] * vis_tri.shape[0], dtype=np.uint32))
                for name, field in data:
                   grid.add_pointdata(DataArray(name, field,
                   vector_format=VF_LIST_OF_COMPONENTS))

                from os.path import exists
                path_to_file = pjoin("vis_files", file_name)
                if exists(path_to_file):
                   raise RuntimeError("output file '%s' already exists"
                   % path_to_file)

                #outf = open(path_to_file, "w")
                #AppendedDataXMLGenerator(compressor)(grid).write(outf)
                #outf.close()

                vis_X_1 = X_1_gpu.get()[b_q:b_q+N_teth]
                vis_X_2 = X_2_gpu.get()[b_q:b_q+N_teth]
                vis_X_3 = np.zeros((vis_X_1.shape))
                vis_X = np.hstack((vis_X_1,vis_X_2,vis_X_3)).copy()
                file_name = "X_teth-%04d.vtu"%t_step
                compressor = None

                grid = UnstructuredGrid(
                     (N_teth, DataArray("points", vis_X, vector_format=VF_LIST_OF_VECTORS)),
                     cells=np.arange(N_teth, dtype=np.uint32),
                     cell_types=np.asarray([VTK_VERTEX] * N_teth, dtype=np.uint8))
                for name, field in data:
                   grid.add_pointdata(DataArray(name, field,
                   vector_format=VF_LIST_OF_COMPONENTS))

                from os.path import exists
                path_to_file = pjoin("vis_files", file_name)
                if exists(path_to_file):
                   raise RuntimeError("output file '%s' already exists"
                   % path_to_file)

                #outf = open(path_to_file, "w")
                #AppendedDataXMLGenerator(compressor)(grid).write(outf)
                #outf.close()

                vis_X_1 = X_1_gpu.get()
                vis_X_2 = X_2_gpu.get()
                #vis_X_1 = X_1_gpu.get()[0:b_q]
                #vis_X_2 = X_2_gpu.get()[0:b_q]
                vis_X_3 = np.zeros((vis_X_1.shape))
                vis_X = np.hstack((vis_X_1,vis_X_2,vis_X_3)).copy()
                file_name = "edges-%04d.vtu"%t_step

                grid = UnstructuredGrid(
                    (b_q+N_teth+b_add,
                    DataArray("points", vis_X,
                    vector_format=VF_LIST_OF_VECTORS)),
                    cells=edges_array.reshape(-1),
                    cell_types=np.asarray([VTK_LINE] * (edges_array.shape[0]), dtype=np.uint32))

                for name, field in point_data:
                    grid.add_pointdata(DataArray(name, field,
                    vector_format=VF_LIST_OF_COMPONENTS))

                for name, field in cell_data:
                    grid.add_celldata(DataArray(name, field,
                    vector_format=VF_LIST_OF_COMPONENTS))

                from os.path import exists
                path_to_file = pjoin("vis_files", file_name)
                if exists(path_to_file):
                   raise RuntimeError("output file '%s' already exists"
                   % path_to_file)

                #outf = open(path_to_file, "w")
                #AppendedDataXMLGenerator(compressor)(grid).write(outf)
                #outf.close()

    queue.finish()
    print 'time spent is', time.time()-start_time, 'seconds'

    f_1_soln = f_1_gpu.get()
    f_2_soln = f_2_gpu.get()
    #f_1_soln2, f_2_soln2 = spreadF(F_1_gpu.get(), F_2_gpu.get(), X_1_gpu.get(), X_2_gpu.get(), b_q, dx, dtheta)
    #print 'err f_1 is', np.max(np.abs(f_1_soln-f_1_soln2))
    #print 'err f_2 is', np.max(np.abs(f_2_soln-f_2_soln2))
    u_1_soln = u_1_gpu.get()
    u_2_soln = u_2_gpu.get()
    X_1_soln = X_1_gpu.get()
    X_2_soln = X_2_gpu.get()
    p_soln = 0.5*(p_gpu.get()+p_last_gpu.get())
    #print u_1_soln
    #print u_2_soln
    #calculate error
    if ii != start_size: 
        errinf1, errinf2, err11, err12, err21, err22, errinfX, err1X, err2X, errinfp, err1p, err2p = err(u_1_soln, u_2_soln, u_1_o, u_2_o, X_1_soln, X_2_soln, X_1_o, X_2_o, p_soln, p_o, f_1_soln, f_2_soln, f_1_o, f_2_o)
        #errinf1, errinf2, err21, err22 = err(u_1_soln, u_2_soln, u_1_o, u_2_o)
        print 'u Errinf1 is', errinf1, '\n u Errinf2 is', errinf2
        print 'u Err11 is', err11, '\n u Err12 is', err12
        print 'u Err21 is', err21, '\n u Err22 is', err22
        print 'X Errinf is', errinfX
        print 'X Err1 is', err1X
        print 'X Err2 is', err2X
        print 'p Errinf is', errinfp
        print 'p Err1 is', err1p
        print 'p Err2 is', err2p
        if ii > start_size+1:
            print 'rinf for u1', -np.log(errinf1/errinf1_o)/np.log(2)
            print 'rinf for u2', -np.log(errinf2/errinf2_o)/np.log(2)
            print 'r1 for u1', -np.log(err11/err11_o)/np.log(2)
            print 'r1 for u2', -np.log(err12/err12_o)/np.log(2)
            print 'r2 for u1', -np.log(err21/err21_o)/np.log(2)
            print 'r2 for u2', -np.log(err22/err22_o)/np.log(2)
            print 'rinf for X', -np.log(errinfX/errinfX_o)/np.log(2)
            print 'r1 for X', -np.log(err1X/err1X_o)/np.log(2)
            print 'r2 for X', -np.log(err2X/err2X_o)/np.log(2)
            print 'rinf for p', -np.log(errinfp/errinfp_o)/np.log(2)
            print 'r1 for p', -np.log(err1p/err1p_o)/np.log(2)
            print 'r2 for p', -np.log(err2p/err2p_o)/np.log(2)
        errinf1_o = errinf1
        errinf2_o = errinf2
        err11_o = err11
        err12_o = err12
        err21_o = err21
        err22_o = err22
        errinfX_o = errinfX
        err1X_o = err1X
        err2X_o = err2X
        errinfp_o = errinfp
        err1p_o = err1p
        err2p_o = err2p
    f_1_o = f_1_soln
    f_2_o = f_2_soln
    u_1_o = u_1_soln
    u_2_o = u_2_soln
    X_1_o = X_1_soln
    X_2_o = X_2_soln
    p_o = p_gpu.get()

#f.close()
